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Abstract 

A  study  of  risk  assessment  for  artillery  fire  randomly  colliding  with 
fixed  wing  aircraft  is  presented.  The  research  lends  itself  to  a  general 
study  of  collision  models.  Current  models  of  object  collisions  fall  under 
one  of  three  categories:  the  historical  model,  the  gas  particle  model,  and 
the  satellite  model.  These  three  vary  in  data  requirements  and  math¬ 
ematical  representation  of  the  impact  event.  The  gas  particle  model  is 
selected  for  its  flexibility  and  robust  estimation.  However,  current  math¬ 
ematical  development  in  the  literature  does  not  include  certain  spatial 
and  dynamic  components  necessary  for  a  general  encounter  (collision) 
model.  These  are  derived  in  this  work.  For  the  specific  application  at 
hand  Quadratic  formulas  estimate  the  ballistic  arc  of  artillery  shells  to 
provide  instantaneous  relative  velocities.  An  extended  Poisson  spatial 
process  is  applied  over  the  relative  volume  within  a  collision  radius  dur¬ 
ing  the  conflict  time  window  to  provide  a  probability  of  collision.  Imple¬ 
mentation  of  the  model  for  military  use  has  been  achieved  via  an  Excel 
spreadsheet  providing  scenario  study  capability  in  real  time.  Results  for 
several  scenarios  are  presented  which  have  been  validated  by  experts. 
These  results  support  current  policy  of  strict  deconfliction. 
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PROBABILISTIC  ESTIMATION  OF  RARE  RANDOM 
COLLISIONS  IN  3  SPACE 


1.  Introduction 

1 . 1  Setting 

The  date  was  March  21,  2001.  The  spring  training  game  between  the  Arizona 
Diamondbacks  and  the  New  York  Giants  was  in  its  7th  inning.  Randy  Johnson 
wound  up  to  pitch  his  famous  hardball,  released,  . . .  and  . . .  the  ball  was  HIT  ...  by 
a  bird!  (Campbell,  2001) 

Random  collisions  occur  all  around  us  every  day.  Meteor  strikes,  particles  of 
air  colliding,  midair  aircraft  collisions,  and  subatomic  particles  smashing  together  in 
nuclear  reactors  are  a  few  examples.  Some  collisions  seem  mundane,  and  others  seem 
like  an  amazing  coincidence.  Some  are  completely  uninteresting,  but  others  grab  our 
attention  with  the  enormity  of  the  event.  Generalizing  on  the  concept  of  collisions, 
there  are  many  interactions  requiring  that  objects  pass  within  some  distance  of  the 
center  of  each  other  to  stimulate  some  event.  Such  events  are  called  “encounters”  - 
where  one  object  encounters  the  other  in  space.  Search  and  rescue  operations  require 
that  the  searchers  pass  within  some  detection  distance  in  order  to  find  the  victims. 
Wireless  networks  require  that  the  client  is  within  a  certain  distance  of  the  hub  in 
order  to  maintain  the  required  signal  strength. 

If  the  motion  of  the  two  objects  can  be  represented  by  some  random  variable, 
then  their  encounter  is  a  random  event  and  can  be  probabilistically  estimated.  Such 
probabilistic  estimates  have  repeatedly  been  requested  regarding  various  encounter 
scenarios  throughout  history.  It  is  argued  that  the  first  question  regarding  random 
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encounters  was  in  the  desire  to  avoid  rainfall  by  running.  However,  that  particular 
scenario  was  not  posed  in  the  literature  until  the  late  20th  century. 

The  earliest  probability  estimate  of  an  encounter  model  dates  back  to  the 
work  of  Rudolph  Clausius  in  the  mid-nineteenth  century.  Clausius  was  advocating 
a  kinetic  theory  of  heat.  At  the  time,  it  was  thought  that  heat  was  some  type  of 
substance,  called  caloric  that  was  transferred  from  one  object  to  another.  Clausius 
believed,  as  is  commonly  held  in  the  scientific  community  today,  that  heat  was  rather 
the  kinetic  energy  of  the  molecules  making  up  the  object.  One  primary  objection  to 
Clausius’  theory  was  the  issue  of  dissipating  vapors:  if  the  gas  particles  really  moved 
as  fast  as  Clausius  claimed,  why  did  it  take  so  long  for  smoke  to  dissipate  from  one 
side  of  the  room  to  the  other.  Clausius’  response  was  to  point  out  that  the  motions 
of  particles  are  not  “free,”  rather  they  are  inhibited  by  collisions  with  other  particles 
randomly  dispersed  throughout  the  space.  He  proceeded  to  demonstrate  that  the 
amount  of  distance  a  particle  can  travel  before  it  randomly  collides  with  another 
particle  is  sufficient  to  prevent  instantaneous  dispersal  of  a  vapor  (Brush,  1965). 

1.2  Problem  Statement 

Understanding  the  issues  regarding  encounter  models  requires  one  to  first  ap¬ 
preciate  the  different  types  of  models  there  are  and  the  variety  of  real-world  encounter 
scenarios  they  represent. 

1.2.1  Background  on  Encounter  models 

Encounter  models  fall  into  one  of  three  categories  discussed  below. 

1.2. 1.1  The  Historical  Model 

The  historical  model  derives  its  name  from  the  general  concept  of  taking  histor¬ 
ical  data  where  similar  incident  objects  have  moved  in  a  space  populated  by  similar 
target  objects  and  considers  the  ratio  of  collision  incidences  to  collision-free  inci- 
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dences  for  an  estimate  of  the  probability  of  future  collisions.  The  model  is  widely 
used  in  many  different  topical  studies.  Mathematically,  it  is  very  simple  and  ex¬ 
tremely  useful  in  cases  where  it  is  impossible  or  undesirable  to  characterize  the  mo¬ 
tion  paths  of  the  objects  whether  by  stochastic  or  deterministic  means.  This  model 
has  been  used  to  estimate  collision  probabilities  in  a  variety  of  scenarios  including 
bird  strikes  on  aircraft  (Bird/Wildlife  Avian  Strike  Hazard  (BASH)  Team,  2008), 
meteor  collisions  on  earth  (Chapman  and  Morrison,  1994),  and  the  general  concept 
was  also  used  in  studies  on  Tomahawk  cruise  missile  fratricides  (Pollack  et  ah,  1997) 
and  satellite  collisions  (Carpenter,  2007). 

In  each  of  the  above  cases  it  would  have  been  extremely  difficult  to  model  the 
motion  paths  of  the  incident  object,  the  target  objects,  or  both.  For  example,  in 
the  case  of  bird  strikes  it  would  be  reasonable  to  expect  the  birds  and  the  pilots  to 
exhibit  avoidance  behavior  when  they  sense  a  collision  is  impending.  On  the  other 
hand,  it  is  much  simpler  to  simply  examine  cases  where  collisions  occurred  and  where 
they  did  not  occur,  and  take  the  ration  of  collision  scenarios  to  all  scenarios  as  the 
probability  of  collision. 

The  historical  model’s  primary  benefits  are  its  simplicity  and  flexibility.  It 
requires  an  input  of  the  number  of  collision  events  and  the  number  of  non-collision 
events.  The  probability  of  collision  is  then  estimated  by 

#  collisions 

Pc  ~  - 

[#  collisions]  +  [#  non-collisions] 

Alternatively,  in  the  case  of  continuous  time  parameters,  it  would  also  suffice  to 
have  several  samples  of  disjoint  time  intervals  tt  for  i  =  1,2, ...  ,n,  each  with  a 
corresponding  count  of  the  collisions  that  occurred  during  that  interval,  N(ti).  Then 
an  estimate  for  the  number  of  events  in  some  future  time  interval,  t* ,  is  given  by 
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multiplying  the  expected  collision  rate  by  the  specified  time  interval: 


Pc 


N{tj)  \ 

tr  nt‘ ) 


r 


The  flexibility  is  demonstrated  by  its  ability  to  deal  with  avoidance  behavior 
as  mentioned  previously,  and  it  is  demonstrated  by  the  model’s  ability  to  accept 
even  quasi- “historical”  data  from  Monte  Carlo  simulations  such  as  the  ones  in  the 
Tomahawk  study  (Pollack  et  al.,  1997).  The  simplicity  allows  for  fast  computation 
of  collision  probabilities  where  the  data  is  readily  available.  For  example,  meteor 
collision  models  consider  the  rate  of  meteor  strikes  and  multiply  this  by  the  prior 
probability  of  a  meteor  being  of  a  certain  size.  This  is  used  in  doomsday  calculations 
for  the  probability  of  a  meteor  destroying  the  human  race  within  the  next  100  years 
(Chapman  and  Morrison,  1994). 

The  largest  difficulty  with  the  historical  model  is  the  standard  question  of  all 
sampled  data:  is  the  data  truly  representative  of  the  whole  population?  One  solution 
for  the  problem  is  to  strip  away  the  quantitative  results  and  simply  give  a  qualitative 
or  comparative  output  (i.e.,  “You  are  more  likely  to  encounter  object  a  than  you  are 
encounter  with  object  b).  This  is  the  approach  taken  by  the  Bird/Wildlife  Aircraft 
Strike  Hazard  Team  in  the  United  States  Air  Force  Avian  Hazard  Advisory  System 
(Bird/Wildlife  Avian  Strike  Hazard  (BASH)  Team,  2008).  Another  answer  to  the 
problem  of  sampling  data  is  by  limiting  the  number  of  scenarios.  The  Tomahawk 
limits  the  study  to  only  the  scenarios  that  are  most  likely  to  result  in  a  fratricide 
and  thus  uses  the  worst-case  scenario  to  provide  an  upper  bound  on  the  probability 
of  collision  for  an  arbitrary  scenario  (Pollack  et  ah,  1997). 

1.2. 1.2  The  Gas  Particle  Model 

The  gas  particle  model  is  a  very  commonly  used  model  for  estimating  collisions 
or  general  encounter  probabilities.  In  this  model  an  incident  particle  moves  through 
the  space  at  a  constant  velocity.  No  information  is  required  or  assumptions  made 
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about  the  location  or  path  of  the  incident  or  target  particle  paths,  except  that  they 
are  uniformly  distributed  throughout  the  space.  All  objects  are  treated  as  point 
particles  and  collisions  occur  when  the  incident  particles  come  within  a  defined  radial 
distance  from  the  target  particles.  Using  the  radial  distance  and  a  constant  forward 
velocity,  the  incident  particles  ”  search”  through  a  measurable  volume  (or  area  for  two 
dimensional  spaces)  in  a  time  interval.  The  target  particle  density  is  then  computed 
over  the  searched  area  to  estimate  the  expected  number  of  collisions.  The  rate  of 
collision  is  the  searched  area  divided  by  the  total  area.  This  rate  is  converted  to  a 
probability  by  using  a  Poisson  arrival  process  for  the  collisions. 

The  model  has  been  heavily  researched  and  applied  to  many  different  scenar¬ 
ios.  As  the  name  suggests,  this  model  is  used  to  examine  gas  particle  interactions 
(Clausius,  1858).  Other  application  domains  for  the  gas  particle  model  include 

•  Sub-atomic  collision  models  (Johnson,  1939) 

•  Minefield  navigation  (Inselmann,  1977;  Kim,  2002) 

•  Orbital  debris  held  calculations  (Jenkin  and  Gick,  2002;  Johnson  et  ah,  2001) 

•  The  question:  “Will  I  stay  drier  by  running  in  the  rain?”  (DeAngclis,  1987) 

•  Bird  strikes  on  windmills  (Ano,  2000) 

•  Aircraft  collision  probabilities  under  “free  flight”  rules  (Alexander,  1970) 

Many  simple  applications  limit  themselves  to  mean  value  analysis  by  assuming 
a  constant  density,  however,  no  such  limitation  is  required  by  the  literature.  Non¬ 
constant  densities  and  constant  density  probabilities  are  handled  by  spatial  point 
processes  (SPPs).  SPPs  are  used  in  a  wide  variety  of  holds  including  forestry,  ge¬ 
ography,  seismology,  cosmology,  spatial  epidemiology,  and  material  sciences  (see  for 
example:  Brix  and  Chadoeuf  (2002);  Comas  and  Mateu  (2007)).  Several  common 
SPPs  are  the  binomial,  Poisson,  Cox,  and  Gibbs  processes.  They  frequently  assume 
stationary  points,  but  not  always  (Mqller  and  Waagepeterson,  2006).  Examples  of 
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non-stationary,  or  “dynamic,”  models  include  epidemiological  models  with  migra¬ 
tion  allowed  (Ball  and  Lynn,  2006)  and  sphere  packing  models  (Stoyan,  2003).  More 
specific  to  the  encounter  modeling  domain,  kinetic  theory  and  orbital  debris  envi¬ 
ronment  models  dynamic  point  processes  of  the  Poisson  variety  to  generate  collision 
probabilities. 

Angelis’  (1987)  solution  to  the  rainfall  question  will  mathematically  map  into 
the  gas  particle  construct,  however  his  particular  solution  incorporates  the  flux  of 
rainfall  against  a  unit  surface.  This  demonstrates  a  key  concept  of  the  model:  relative 
motion.  In  Angelis’  formulation,  he  takes  the  frame  of  reference  of  a  stationary 
observer  and  the  particles  are  oncoming,  whereas  many  other  applications  use  the 
frame  of  reference  where  the  observer  is  moving  in  a  field  of  stationary  particles.  In 
many  cases  it  is  conceptually  easier  to  work  in  the  latter  frame  of  reference,  which 
is  the  approach  we  will  use.  Note  that  relative  velocities  are  reflexive,  so  the  two 
formulations  are  equivalent. 

The  strength  of  the  gas  particle  model  is  its  wide  applicability  in  many  domains. 
With  relatively  simple  calculations,  this  model  has  been  adjusted  for  many  different 
assumptions  on  the  locations  and  paths  of  incident  and  target  particles,  allowing  for 
relatively  flexible  applications.  More  advanced  adjustments  to  the  basic  model  can 
be  used  to  account  for  target  particle  motions,  number  of  dimensions,  and  relativistic 
matter/energy  considerations  in  the  case  of  sub-atomic  physics. 

The  gas  particle  model’s  weaknesses  for  representing  encounter  scenarios  is 
that  it  can  become  difficult  to  model  dynamic  situations.  Furthermore,  the  model 
requires  that  any  object’s  motion  is  independent  of  the  motions  of  other  bodies  in 
the  space.  No  targeting  or  avoiding  behavior  is  permitted. 

1.2. 1.3  The  Satellite  Collision  Model 

This  third  model  requires  more  information  than  either  of  the  previous  two. 
Its  primary  field  of  application  is  in  estimating  the  probability  that  two  satellites 
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in  orbit  will  collide  (Alfano,  2006;  Alfriend  et  al.,  1999;  Carlton-Wippern,  1993; 
Chan,  1997).  Another  common  application  is  in  airspace  deconfliction  by  air  traffic 
controllers  (Chiang  et  ah,  1997).  In  both  of  these  cases,  there  is  complete  data  for 
estimating  positions,  velocities,  and  pairwise  covariances  of  the  object  motion  paths. 
This  model  would  not  work  well  for  predicting  interactions  where  the  exact  motion 
paths  of  the  particles  are  erratic  and  generally  unknown;  therefore,  it  has  not  been 
applied  to  as  wide  a  variety  of  domains  as  the  other  two  because  of  its  enormous 
data  requirements. 

The  collision  model  requires  full  information  on  the  perceived  location  and  ex¬ 
pected  path  of  every  object  under  consideration.  It  then  uses  geometric  techniques 
to  compute  the  points  of  closest  approach  for  every  incident-target  object  pair.  Lay¬ 
ering  a  Gaussian  error  distribution  at  the  point  of  closest  approach,  an  upper  bound 
on  the  probability  of  collision  is  estimated  by 


where  V  is  the  volume  of  space  within  which  the  encounter  is  expected  to  occur,  C  is 
the  combined  covariance  of  the  pair  of  objects  and  f  is  the  collision  distance  (Chan, 
2004). 

1.2.2  Discussion 

Between  these  three  modeling  techniques  a  large  gap  exists  for  encounter  sce¬ 
narios  which  have  data  limitations,  where  encounters  are  very  rare,  and  where  the 
physical  behavior  is  complex.  In  this  research  we  will  extend  the  gas  particle  model 
with  a  view  towards  filling  in  this  gap.  Specifically,  some  examples  of  scenarios  that 
break  the  assumptions  of  existing  gas-particle  encounter  models  include: 

•  Speeds  of  objects  may  be  dependent  on  the  objects’  locations  in  space 
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•  Objects  are  not  uniformly  distributed  throughout  the  space 

•  The  number  of  encounter-able  objects  moving  in  the  space  is  bounded  by  some 
positive  integer  n 

•  Objects  may  enter  or  leave  the  collision  space 

1.3  Objective 

The  goal  of  this  research  is  to  advance  encounter  model  probability  theory  to 
provide  collision  risk  estimates  for  more  general  spatial  and  dynamic  cases.  The  sub 
objectives  to  be  completed  are: 

1.  Research  modern  point  process  theory  for  feasible  enhancements  to  extend  the 
encounter  model  theory. 

2.  Apply  the  enhanced  model  to  a  problem  scenario  to  provide  a  proof  of  concept 
for  the  model. 

1 . 4  Scope 

The  model  developed  in  this  research  could  apply  to  a  number  of  different 
scenarios.  The  application  selected  for  this  research  is  one  of  great  military  interest: 
a  probability  estimation  that  artillery  fire  may  randomly  encounter  aircraft  in  midair. 
Figure  1.1  provides  the  problem  scenario.  The  artillery  unit  is  firing  indirect  fires  (i.e. , 
the  gunner  cannot  see  the  target)  that  arc  through  an  altitude  range  where  aircraft 
are  flying.  Clearly  such  a  scenario  presents  risk  to  the  aircraft,  but  Musselman  (2008) 
explains  further  that  if  pilots  are  concerned  about  being  hit  by  friendly  artillery  fire, 
they  may  not  be  as  careful  with  their  munitions  as  they  would  otherwise  be,  since 
they  may  be  distracted  by  the  extra  threat  to  their  aircraft.  The  threat  to  aircraft 
from  the  ground  artillery  increases  risks  for  both  the  aircrews  and  the  ground  forces. 
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Figure  1  -  Problem  scenario 


Figure  1.1  The  Problem  Scenario 

Current  doctrine  prevents  this  from  ever  occurring  by  requiring  commanders 
to  deconflict  the  airspace.  Deconfliction  consists  of  using  fire  support  coordinating 
measures  (FSCMs)  to  allocate  portions  of  the  airspace  for  certain  use  during  specified 
time  intervals.  Common  FSCMs  include  lateral  separation  (i.e.,  you  stay  on  your 
side  of  this  line,  and  I’ll  stay  on  my  side),  altitude  separation  (e.g.  the  aircraft  flies 
above  the  artillery  fires),  or  time  separation  (JP3-09.3,  2003).  The  kill  box  structure 
is  an  example  of  a  more  complex  FSCM  that  is  currently  used  in  military  operations 
(JP3-09.34,  2005). 

Despite  the  fact  that  joint  doctrine  requires  the  application  of  FSCMs,  risk 
analysis  for  artillery  fires  colliding  with  aircraft  in  a  conflicted  airspace  is  increas¬ 
ingly  becoming  a  point  of  interest  for  mission  planners,  decision  makers,  and  com¬ 
manders.  In  December  2007,  the  Center  for  Army  Lessons  Learned  published  a 
requirement  for  an  analytical  model  for  assessing  the  risk  of  aircraft  fratricide  due  to 
fires  (Neuenswander,  2008).  Bethea  and  Herbranson  (2008)  point  out  that  several 
recent  exercises  on  the  Korean  peninsula  highlight  the  difficulty  of  deconfliction  in 
small  airspaces  that  are  heavily  congested.  Deconfliction  in  the  Korean  peninsula 
places  limits  on  the  commander’s  options,  ability  to  mass  his  forces  effectively,  and 
conduct  timely  responses  to  new  developments  on  the  battlefield.  Smith  (2006)  also 
identifies  that  in  aerial  interdiction  operations  the  FSCMs  can  provide  sanctuary  to 
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the  enemy,  prevent  the  prosecution  of  time  sensitive  targets,  and  limit  the  Comman¬ 
der’s  and  pilot’s  tactical  options. 

Under  these  limitations,  the  commander  is  left  in  a  bind  where  he  must  choose 
between  one  of  two  options: 

1.  Accept  no  risk  to  aircraft  from  friendly  surface  fires  by  deconfliction.  This  will 
result  in  limitations  as  described  above. 

2.  Accept  an  unknown  risk  that  his  aircraft  will  be  hit  by  friendly  surface  fire  by 
choosing  not  to  deconflict  the  airspace.  This  will  provide  more  flexibility  since 
many  of  the  limitations  described  above  are  not  applicable. 

Clearly,  an  unknown  risk  makes  the  decision  very  difficult  to  commanders. 

There  have  been  two  prior  efforts  at  determining  this  unknown  risk.  The 
first  was  conducted  by  the  National  Defense  Industrial  Association  (NDIA)  in  2001 
(Strike,  Land  Attack  and  Air  Defense  (SLAAD)  Committee  of  the  National  Defense 
Industrial  Association  (NDIA),  2001).  Their  study  was  conducted  on  behalf  of  the 
Office  of  the  Chief  of  Naval  Operations.  The  study  covered  many  aspects  of  Joint 
Fires  including  various  methods,  tools,  and  issues  with  the  practice.  As  a  sub¬ 
objective,  the  NDIA  team  conducted  a  simulation  to  test  the  validity  of  the  “Big 
Sky  -  Little  Missile”  theory  that  fratricide  risks  are  negligibly  small.  A  number  of 
important  details  are  omitted  from  the  study  report,  including  confidence  intervals 
and  numbers  of  replications  of  the  simulation.  However,  they  conclude  that  the  risk 
that  aircraft  will  be  hit  by  friendly  fire  is  1%,  which  is  not  a  negligible  risk.  No 
sensitivity  analysis  seems  to  have  been  performed  on  the  results  to  test  different 
scenarios. 

In  2008  the  Air  Force  Doctrine  Development  and  Education  Center  (Bethea 
and  Herbranson)  started  a  probabilistic  model  that  reduces  the  airspace  to  the  ver¬ 
tical  plane  that  cuts  through  the  gun-to-target  line  (GTL).  The  plane  is  then  par¬ 
titioned  into  aircraft-sized  rectangles.  The  ballistic  path  is  approximated  by  an 


1-10 


isosceles  triangle  on  the  plane.  The  risk  estimate  is  given  by  the  ratio  of  rectangles 
intersected  by  the  triangle  to  the  total  number  of  rectangles.  This  is  a  much  cruder 
approximation  than  we  propose  to  provide  with  our  research. 


To  date,  no  published  model  exists  for  probabilistic  estimation  of  the  dynamic 
collision  likelihood  of  artillery  fires  on  fixed  wing  aircraft.  Specifically,  a  model  is 
needed  that  can  account  for 

•  Uncertainty  in  the  artillery  location 

•  Uncertainty  in  the  aircraft  location/heading 

•  Uncertainty  in  the  target  location 

•  Uncertainty  in  the  shell  speed 

•  Uncertainty  in  the  aircraft  speed 

•  Collisions  from  any  possible  angle  (head-on,  from  behind,  from  above,  etc.) 

•  Duration  of  airspace  confliction 

•  Number  of  aircraft  in  the  airspace 

•  Aircraft  altitude 

•  Artillery  max  ord  (highest  height  achieved  by  the  fired  shell) 

•  Size  of  the  operational  area  (surface  area  from  a  top-down  view) 

1.5  Overview 

Chapter  2  provides  a  theoretical  development  of  an  enhanced  gas  particle  model 
where  the  object  speeds  are  spatially  determined  and  their  spatial  distribution  fol¬ 
lows  a  generalization  of  the  spatial  Poisson  point  process.  Chapter  3  presents  an 
application  of  the  enhanced  model  to  the  proposed  study  for  aircraft  fratricides, 
along  with  sensitivity  analysis  and  numerical  results.  Chapter  4  concludes  the  thesis 
by  highlighting  the  key  research  hirelings  and  some  areas  for  future  research. 
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2.  Theoretical  Development 

Advancing  Encounter  Models  with 
Modern  Point  Process  Theory 


2. 1  Introduction 

We  consider  the  general  class  of  encounter  models:  representations  of  objects 
moving  in  space  as  they  randomly  encounter  other  objects  in  the  space.  It  is  assumed 
that  the  space  is  M3.  Common  applications  of  encounter  modeling  techniques  include 
kinetic  theory  (i.e. ,  motions  and  collisions  of  gas  particles)  (Brush,  1965),  bird  strikes 
(Bird/Wildlife  Avian  Strike  Hazard  (BASH)  Team,  2008),  airspace  deconfliction 
(Alexander,  1970;  Pollack  et  ah,  1997),  and  satellite  collision  models  (Chan,  2004). 
The  result  is  a  more  general  spatial  point  process  (SPP)  which  is  fully  dynamic. 

2.1.1  The  Difficulty  on  Dynamic  SPPs 

Mqller  and  Waagepeterson  (2004)  define  an  SPP  X  on  a  space  S  C  Mn  as  a 
measurable  mapping  on  some  probability  space  (Q,  IF,  /)  to  ( Nif ,  where  Nif 

is  the  set  of  locally  finite  point  configurations  in  S  and  <t(A/)j)  is  it’s  o- algebra.  The 
distribution  of  X  follows  Px{F)  =  P{{u  £  ^  :  X(co)  G  F})  for  any  F  G  Nif.  They 
go  on  to  point  out  that  every  SPP  can  be  uniquely  identified  by  its  corresponding 
counting  process  on  the  collection  of  closed  Borel  sets  in  S,  which  is  the  approach 
taken  here.  We  define  a  random  locally  finite  point  configuration  of  N  points  on  a 
set  S  C  M3  with  the  state  space 

B  =  {Bi,  B2,...,Bn:  Bi  G  S} 


2-1 


Then  we  say  that  a  SPP  is  binomial  if  for  any  W  C  S,  it  holds 


Nx{W)  =  \W  n  B\  ~  Binomial ( | B |,  f(W)) 


for  some  density  function  /  on  S.  Similarly,  the  Poisson  SPP  is  defined  such  that 
for  any  W  C  S,  \B  D  W\  ~  Poisson(A .fJ>(W))  and  conditional  on  \B  D  W\  =  N,  the 
process  inside  W  is  binomial. 

In  most  applications  of  SPPs,  the  points  are  assumed  to  be  stationary  (Mpller 
and  Waagepeterson,  2006),  which  may  not  be  a  problem  for  certain  encounter  models 
(see  Kim  (2002)  for  example).  However,  in  many  encounter  scenarios  all  the  particles 
are  in  motion,  rendering  the  system  dynamic. 

The  common  technique  to  solve  such  systems  is  to  change  the  frame  of  reference 
so  that  only  one  object  is  in  motion  and  all  the  rest  are  stationary.  Loeb  (1934) 
provides  a  proof  that  the  expected  relative  velocity,  vrei  between  two  objects,  A  and 
B,  moving  in  M3,  provided  that  the  directions  of  motion  for  A  and  B  are  independent 
and  at  least  one  is  uniformly  distributed,  is 


E  e[vrei(vA,vB,0)] 


VA  +  Sj,  VA>VB 
VB  +  i VA<VB 


(2.1) 


where  vA  and  vB  are  the  speeds  of  A  and  B  respectively,  and  9  is  the  angle  between 
their  motion  vectors.  As  a  result,  the  average  rate  of  encounters  in  increased.  Un¬ 
fortunately,  since  this  is  a  rate  increase  it  must  be  evaluated  at  infinitesimal  level. 
Comas  and  Mateu  in  their  simulation  of  a  Gibbs  process  highlight  the  difficulty  of 
examining  such  dynamical  systems  using  discrete  time-step  methods  (Comas  and 
Mateu,  2007). 

While  the  Poisson  process  naturally  scales  to  account  for  changes  at  the  in¬ 
finitesimal  level,  non-Poisson  distributions  such  as  the  binomial  do  not  scale  so  easily. 
Evaluating  instantaneous  increases  in  interaction  probabilities  in  non-Poisson  spa- 
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tial  distributions  is  a  non-trivial  exercise.  In  fact,  if  the  scaling  is  attempted  in  a 
naive  method,  it  can  result  in  probabilities  greater  than  unity,  as  will  be  seen  later. 
Therefore,  while  quite  a  few  dynamic  models  exist  that  employ  the  Poisson  pro¬ 
cess  or  its  derivatives  (specifically  the  Gibbs  process)  to  examine  dynamic  encounter 
systems  (Baddeley  et  ah,  2006;  Clausius,  1858),  to  our  knowledge  no  studies  have 
been  performed  to  examine  dynamic  SPP  systems  with  non-Poisson-based  counting 
processes.  Population  epidemic  models  with  migration  might  be  considered  an  ex¬ 
ception;  however  those  cases  consider  migration  on  a  lattice  structure  rather  than  in 
a  compact  space  (Ball  and  Lynn,  2006;  Hagenaars  et  al.,  2004). 

2.1.2  The  Extended  Poisson  Process 

A  natural  solution  to  this  problem  of  infinitesimal  estimates  on  non-Poisson 
distributions  is  provided  by  the  extended  Poisson  process,  adapted  from  the  chrono¬ 
logical  M  domain  to  the  spatial  M3  domain.  The  core  concept  of  the  extended  Poisson 
process  is  to  use  Markov  birth-death  processes  to  model  counting  processes  that  are 
over-  or  under-dispersed  with  respect  to  the  Poisson  distribution.  The  extended 
Poisson  process  uniquely  features  transition  rates  that  are  functionally  dependent 
on  the  current  state.  For  example,  the  well  known  pure  birth  Yule  process  is  an 
extended  Poisson  process  with  a  linearly  increasing  arrival  rate. 

The  extended  Poisson  process  was  developed  in  the  late  1980’s  through  the 
early  1990’s  in  a  series  of  papers  primarily  between  M.  Faddy,  F.  Ball,  and  P.  Don¬ 
nelly  (Ball,  1995;  Ball  and  Donnelly,  1987;  Brown  and  Donnelly,  1993;  Faddy,  1985, 
1990,  1994).  The  culmination  of  this  chain  of  research  was  when  Faddy  showed  that 
any  discrete  distribution  with  finite  support  can  be  represented  by  a  Markov  death 
process  where  the  transition  rates  depend  on  the  current  state  (1998). 

Since  the  discovery  of  the  extended  Poisson  process,  further  refinement  has 
been  in  the  realm  of  advancing  statistical  data  analysis  using  the  extended  Poisson 
process.  Parameter  estimation  techniques  for  the  extended  Poisson  process  model 


2-3 


include  maximum  log-likelihood  calculations  for  the  Q  matrix  (Podlich  et  ah,  1999), 
iteratively  re-weighted  least  squares  algorithms  (Toscas  and  Faddy,  2005)  and  non- 
parametric  sequence  algorithms  (Podlich  et  ah,  2004).  It  is  interesting  to  note  that 
Quigley  and  Walls  solve  a  similar  problem  where  the  arrival  rate  is  both  dependent 
on  the  state  of  the  system  (number  of  arrivals)  and  the  current  time  by  using  a 
bootstrapping  technique  (2005).  In  application  of  statistical  analysis,  the  extended 
Poisson  process  has  been  heavily  used  in  studies  of  biological  processes  (Faddy,  1995, 
1997;  Faddy  and  Smith,  2008,  2005;  Faddy  et  ah,  2001).  We  are  aware  of  no  major 
further  developments  of  the  extended  Poisson  process  as  a  feature  of  probability 
theory. 

With  the  bulk  of  research  in  extended  Poisson  processes  going  into  parame¬ 
ter  estimation,  the  application  of  extended  Poisson  processes  to  spatial  processes 
is  largely  untouched.  Faddy  performed  a  short  study  of  spatial  data,  but  overtly 
bypassed  a  full  development  of  the  formulaic  application  of  the  extended  Poisson 
process  to  the  spatial  parameters  of  the  problem  (Faddy  and  Smith,  2008).  We 
find,  however,  that  using  the  Lebesque  measure  and  integration  allows  us  to  easily 
port  the  formulas  from  the  time  domain  to  the  spatial  domain.  Most  of  the  key 
results  carry  over  without  any  fundamental  change.  However,  in  higher-dimensional 
domains,  we  find  a  new  usefulness  of  the  extended  Poisson  process:  the  ability  to 
smoothly  handle  increases  in  rates  due  to  the  relative  frames  of  reference  used  in 
encounter  models. 

2.1.3  Outline 

Section  2  lays  the  groundwork  for  our  discussion  by  defining  the  necessary 
terminology  and  assumptions.  Section  3  introduces  the  concept  of  motion  to  develop 
a  dynamic  SPP.  The  section  concludes  with  the  description  of  an  intuitive  step 
which  would  seem  to  allow  for  estimates  of  encounter  probabilities,  but  which  in 
fact  leads  to  a  paradoxical  relationship.  The  solution  is  presented  in  Section  4  with 
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the  introduction  of  the  Extended  Spatial  Poisson  Point  Process  and  an  example  of 
the  model  is  presented  for  a  dynamic  SPP  scenario.  Some  concluding  thoughts  are 
provided  in  Section  5. 

2.2  Definitions,  Assumptions,  and  Scenario  Description 

We  consider  the  probability  that  a  point  A  will  encounter  some  point  B  from 
a  collection  B.  Mathematically,  the  number  of  encounters  is  given  by  a  counting 
process,  Nx,  on  the  r-ball  of  A,  f3r(A)  —  (s  G  S  :  |s-—  A\  <  r}  (using  the  standard 
1 2-norm). 

Nx  =  ^2lMA)(B) 

BeB 

where  I  is  the  indicator  function  of  whether  B  is  in  (3r{A).  Since  this  counting 
process  is  involved,  we  introduce  the  following  assumption  that  allows  us  to  apply 
SPP  theory. 

Assumption  2.1  X  is  a  binomial  SPP  with  uniform  density  function.  Thus,  if  B 
is  distributed  according  to  X,  \B  fl  W\  ~  Binomial(N,  p(W) / p(S)) ,  where  N  =  \B\. 

We  use  the  standard  Lebesgue  measure  p  which  is  assumed  to  be  in  3-space.  Later 
in  the  development  we  will  require  Lebesgue  measures  in  lower  dimensions,  at  which 
time  we  will  subscript  the  measure  to  indicate  the  number  of  dimensions. 

2.2.1  Defining  locations  and  motion 

The  core  of  encounter  models  is  the  concept  of  motion.  Motion  is  commonly 
known  to  be  the  change  in  position  over  time.  Frequently,  the  velocity  of  some  object 
X  is  denoted  by  a  vector,  vx  =  (v,  6,  0),  where  v  G  [0,  oo),  9  G  [0,  27r],  and  0  G  [0,  n\. 
Of  these,  the  most  important  in  this  paper  is  the  scalar  v  which  we  define  as  follows: 
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Definition  2.1  The  speed  of  an  object  in  S  is  given  by  0  a  mapping  0  :  [0,  oo)  — > 
[0,  cx))  given  by  some  constant  v.  The  domain  is  some  time  duration  t,  and  the  image 
is  the  measure  of  some  distance,  d  =  fi\ [d0,df\.  Thus  t>(£)  =  vt  =  d  =  /^([d  —  0,  df\). 

The  above  function  definition  may  be  read  “d  is  the  distance  traveled  by  moving 
at  speed  v  for  time  f.”  Note  that  the  restriction  to  the  non-negative  reals  eliminates 
consideration  of  negative  times  or  negative  speeds.  Throughout  this  article  we  will 
define  v  as  a  function  space,  s  G  S.  Thus,  the  reader  may  interpret  v(s  e  S)  as 
d  =  v(s,  t )  =  v(s)t,  which  is  the  distance  that  would  be  traveled  by  moving  at  speed 
v(s)  for  time  t. 

In  the  dynamic  encounter  problem  considered  here,  both  A  and  the  elements 
of  B  can  be  in  motion  (i.e.  have  positive  speed).  Most  of  kinetic  theory  literature 
(hard-body  particle  collisions  on  Poisson  SPPs)  examines  cases  where  the  elements 
of  B  have  speeds  that  are  independent  and  identically  distributed  (see  for  example 
Clausius  (1858);  Maxwell  (1858)).  To  the  authors’  knowledge,  no  prior  work  exists  in 
the  literature  that  examines  encounter  cases  where  the  speed  of  an  arbitrary  element 
of  B  is  a  function  of  the  element’s  position  in  S.  Examples  of  this  phenomenon 
abound  in  the  natural  environment:  traffic  lanes  have  different  speed  distributions, 
ballistic  projectiles  vary  their  velocity  by  altitude,  aircraft  at  different  altitudes  fly 
at  different  speeds,  and  satellite  speeds  change  between  apogee  and  perigee  of  orbit. 
Note  that  in  some  of  these  cases  the  velocities  are  deterministic  and  in  others  they 
are  stochastic  in  nature. 

It  would  be  possible  to  examine  cases  where  the  velocities  are  random  functions 
of  the  environment  using  random  fields,  but  we  believe  the  core  concepts  can  be 
adequately  expressed  in  deterministic  form.  Thus,  we  consider  the  specific  case  where 
the  function  for  speed  is  deterministic  and  the  direction  is  uniformly  distributed  in 
M3. 
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Assumption  2.2  The  velocities  of  all  the  elements  B  e  B  are  i.i.d.  with  Vg  = 
(r?g(s),  9,  (f)  where  vb(s)  :  S  — >  [0,  oo)  is  deterministic,  9  ~  Uniform^),  2n],  and 
(f)  ~  Uniform[0,7r]  (where  ~  means  “distributed  as”). 

Note  2.1  We  recognize  that  Assumptions  2.2  and  2.1  can  place  certain  limitations 
on  vb{s )  and  X{t).  For  example,  the  following  scenario  presents  a  contradiction: 

{1  s  =  s* 

Vt  e  At 

0,  s^s* 

since  the  system  initializes  with  all  elements  of  B  located  at  s* .  As  the  system  evolves, 
s*  is  transient,  and  every  other  location  is  absorbing,  which  contradicts  that  B  e  B 
is  located  at  s*  with  probability  1.  Thus,  by  the  statement  of  these  two  assumptions, 
it  is  further  implicitly  assumed  that  no  such  contradictions  are  possible. 


Lemma  2.1  If  B  is  distributed  according  to  X  in  S,  then  the  probability  that  B  is 
in  some  region  W  C  S  is  given  by  f(W)  =  E w  [/(s)]  /^(W),  where  E^[/]  means  the 
expectation  of  /  over  W . 

Proof.  The  result  follows  trivially  from  the  definition  of  the  expectation: 


mn 

A*(W0 


Having  addressed  the  motion  of  the  elements  of  B,  we  now  consider  the  motion 
of  A,  which  is  intimately  related  to  the  encounter  region.  Since  encounters  are 
spatially  defined,  the  occurrence  of  an  encounter  depends  on  the  location  of  the 
encounter ing  and  the  encountered.  In  this  problem  formulation,  A  is  the  encountering 
object,  and  the  elements  of  B  are  the  encountered  objects.  The  locations  and  motions 
of  the  encountered  (or  unencountered)  elements  of  B  are  given  by  Assumptions  2.1 
and  2.2  respectively.  We  track  A’s  location  by  its  path  as  follows: 


2-7 


Definition  2.2  A  path,  p,  is  a  continuous  measurable  mapping  from  [to,  tn]  C  Mq 
into  S  where  Jp(p(t))=  vA(t). 

Here  Jp(p(t))  is  the  3-dimensional  Jacobian  of  the  path  p  from  A’s  location  at  time 
t.  Definition  2.2  allows  for  either  stochastic  or  deterministic  motion  of  A.  The 
requirement  on  the  Jacobian  limits  the  path  length  to  what  would  be  achievable 
with  the  given  velocities: 


/  vA(t)dt  =  /  p1(vA(t))dt=  /  vA(t)dt  (2.2) 

Jt0  J  Jto  Jt0 

Clearly,  in  order  to  define  a  path  pA(t)  according  to  Definition  2.2,  we  require 
the  velocity  distributions.  Previous  work  assumes  that  A’s  motion  is  equal  in  dis¬ 
tribution  to  that  of  B.  While  this  is  certainly  allowable,  this  is  not  required  by  the 
mathematical  formulation  of  the  encounter  model.  In  fact,  this  model  can  work  for 
A  with  any  distribution  on  vA  (including  that  of  no  motion  at  all).  For  simplicity 
we  assume  that  the  speed  of  A  is  i.i.d.  at  all  moments  in  time.  This  means  from 
Equation  (2.2)  that  E[/xi  (p)  =  E[uJ4](tn  —  t0).  We  make  no  assumptions  on  the 

directional  components  of  vA. 

Assumption  2.3  vA(t)  =  (vA,6(t),(f)(t))  where  Oft)  ~  fg(t)  and  <f  ~  for  some 
random  variables  6  and  <f>  mapping  t  to  the  probability  spaces  ([0,  27t],  ^[ogTr],  fe(t)) 
and  ([0, 7r],  ^-’[0)7r], /^(t))  respectively  (where  is  the  Borel  a-algebra  on  Q ). 

2.2.2  Defining  an  encounter  region 

Having  located  A  in  space,  we  can  now  proceed  to  define  the  region  in  which 
the  elements  of  B  are  to  be  encountered.  We  do  this  using  a  closed  ball: 

Definition  2.3  Define  /3(d,r)(s)  as  the  closed  r-ball  of  d  dimensions  (d  G  {0, 1,  2,  3} ) 
at  s  G  S  by  /3^rfis)  =  (s'  G  S  :  p(s  —  s')  <  r}  for  some  r  G  M+ 
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Similarly  to  the  definition  of  /j,  we  will  forego  the  use  of  dimensional  subscripts 
for  3  dimensions.  Since  SCI3  we  know  that  j3r{s)  is  unique.  Similarly,  /3(o,r)(s)  =  s 
is  unique.  For  d  =  1  or  d  =  2,  /3^,r)(s)  is  not  unique  in  M3.  Thus,  when  we 
use  /3(2,r)  (s)  below,  we  place  the  further  constraint  on  the  2-dimensional  ball  that 
P(2 ,r)(p(t))  -L  Jp(p(t))i  which  makes  it  unique  in  M3.  For  defining  the  encounter 
region  in  M3,  the  ball  definition  works  for  a  single  point,  but  since  A  may  be  moving 
in  S,  we  must  extend  this  definition  to  include  all  points  covered  by  A  during  its 
sojourn  along  p(t). 

Definition  2.4  Let  A  move  along  pit)  for  t  G  At  =  [to,tn].  Then  let  Wr(p,  At)— 
{s  G  S  :  3t  G  At  =  [t0,t„]  3  s  G  f3(d,r)(p(t)) ■  We  call  Wr(p,  At)  the  encounter 

window  of  A  along  pit)  for  t  G  At. 

Using  the  window  construct  is  critical  for  encounter  models.  From  Lemma 
2.1,  using  p(Wr(p,  At),  we  can  find  the  density  function  integrated  over  W,  which  is 
used  by  the  binomial  process  X  (given  in  Assumption  2.1)  to  compute  the  random 
counter  of  the  number  of  elements  in  B  encountered  by  A  during  At. 

In  order  to  find  W,  we  break  with  the  standard  convention  in  kinetic  theory 
of  using  the  disks  perpendicular  to  Va,  and  integrate  over  the  path.  Instead,  we 
will  use  a  “leading  edge”  construct.  The  result  is  the  same  “tube”-shaped  encounter 
window  of  measure  pd( W)  =  Af(<i-i)(/3(d-i,r))Afi(p)  under  certain  assumptions  (see 
Assumption  2.4  below)  on  W  that  prevent  it  from  doubling  back  on  itself  with  a 
positive  Us  measure.  We  begin  by  defining  the  leading  edge.  The  reason  for  the  use 
of  the  leading  edge  as  opposed  to  the  disk  will  be  explained  at  the  end  of  this  section. 

Definition  2.5  Let  A  move  along  p  in  S  over  some  time  At  =  [to,tn].  At  any  time 
t  G  At,  let  V(s)  denote  the  plane  perpendicular  to  VA(t)  passing  through  some  point 
s  in  the  image  of  p.  We  define  the  leading  edge,  LEp(s),  as  the  intersection  of  (3=(s) 
with  the  closed  half- space  bounded  by  V(s)  and  on  the  same  side  of  Vis)  as  VA(t), 
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AT) 

Figure  2.1  The  leading  edge  in  2  dimensions 

where  /3Tdr^(s)  is  the  boundary  of  the  r-ball  given  by  an  equality  constraint  on  the  ball 
definition. 

Figure  2.1  gives  a  2-dimensional  represetation  of  the  concept  of  the  leading 
edge.  We  constrain  the  path  to  prevent  the  encounter  window  from  doubling  back 
on  itself.  Clearly,  this  cannot  be  done  with  balls,  since  the  balls  for  two  points 
{t,  t  +  dt}  intersect  each  other  at  uncountably  many  points.  However,  we  can  achieve 
the  desired  constraint  using  the  leading  edge. 

Assumption  2.4  If  A  moves  according  to  path  p  C  S  over  time  At  and  the  en¬ 
counter  radius  is  some  r  >  0,  then  for  any  ta,tb  C  At  where  ta  h  it  holds  that 

LEp(p(ta))  fl  Pr(p(tb))  =  0. 

Lemma  2.2  p  is  one-to-one. 

Proof.  Assume  the  corollary  is  false.  Then  3 ta,  tb  G  At  where  ta  and 
p(ta)  =  pitfi).  Clearly,  this  would  mean  LEp(p(ta ))  C  (3r(p(ta))  =  /3r(p(tb)),  which 
contradicts  Assumption  2.4.  ■ 

Note  2.2  It  can  be  shown  that  Assumption  2.f  holds  if  p  has  curvature  less  than  -. 

Lemma  2.3  Suppose  \B\  =  N  and  the  elements  B  are  stationary  and  distributed 
in  S  according  to  the  binomial  SPP  X .  Then  let  A  move  along  p  for  some  time 
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At.  Suppose  further  that  it  is  known  that  for  some  connected  interval,  8\ t  C  At 
with  t0  G  8t\ ,  that  n  <  N  elements  of  B  were  within  radius  r  from  A  at  some 
time  during  S\ .  Then  the  probability  of  encountering  further  elements  in  the  time 
interval  S2t  =  At\S\t  is  given  by  the  formula:  P[Nx(Wr(p,  At))\N(Wr(p,  Sit)  =  n 
and  N(S)  =  N]  ~  ( Binomial(N  —  n ,  f(Wf(p ,  S2 1)))  +  n). 

Proof.  Since  X  is  binomial  and  it  is  known  that  N  —  n  elements  are 
contained  in  S\Wr(p,  5\t),  the  result  follows  from  the  definition  of  X.  ■ 

Next  we  consider  the  new  encounter  region  swept  by  some  particle  A  moving 
along  path  p  during  a  connected  time  interval  5t  C  At.  This  region  is  given  by 
W*(p,8t)—  Wr(p,  8t)\/3r(p(ti))  where  t\  =  inf(ht)  and  if  t\  e  St  we  exclude  the  open 
r-ball,  otherwise  we  exclude  the  closed  ball.  Using  the  leading  edge  this  region  is 
given  by  the  union 

W:(p,St)  =  \jLEp(p(t))  (2.3) 

test 

Since  this  union  is  incremented  as  t  increases,  we  refer  to  it  as  the  incremental 
encounter  window. 

Our  final  lemma  in  this  section  comes  in  light  of  Lemmas  2.1  and  2.3.  With 
the  end  goal  of  finding  Ng(W*(p,5t)),  from  Lemma  2.3  we  are  required  to  find 
f(W*(p,  S2t))  which  with  the  help  of  Equation  (2.3)  is  shown  to  be 

f(w;(p,S2t))  =  f([  LEp(p(t))dp(t)). 

J  62 1 

This  may  still  be  difficult  to  integrate,  but  from  Lemma  2.1  we  can  change  the  right 
hand  side  to  the  product  of  E[f(Wf(p,At\S\t))]  and  p,(JS2tLEp(p(t))dp(t)).  This 
places  us  in  the  position  for  the  final  result  of  this  section: 

Theorem  2.1  Suppose  A  is  moving  along  p  for  some  time  St.  Then  p(Wf(p,  St))  = 

7Tr2VASt. 
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Proof.  We  know  from  Equation  (2.3)  that  p(Wf(p,  St))  =  p(fSt  LEp(p(t))dp(t). 
Then  from  the  properties  of  the  Lebesgue  measure  p 

M  f  LEp(p(t))dp(t))  =  [  p(LEp(p(t))dp(t)) 

JSt  Jst. 

=  f  p2(LE*(p(t)))p1(dp(t)) 

Jst 

=  /  nr2 pi(vA(t)dt)dt 

Jst 

=  7T-r2  /  VAPxidt) 

Jst 

=  nr2vApi(St)  (2.4) 

where  p-2 (LE*(p(t)))  is  the  surface  area  of  LEp(p(t))  projected  onto  the  plane  per¬ 
pendicular  to  va(J),  which  is  a  disk  of  radius  r.  The  projection  is  necessary  because 
the  Lebesque  measure  can  only  be  separated  into  constituent  dimensions  if  those 
dimensions  are  orthogonal.  ■ 

It  is  interesting  to  note  that  Theorem  2.1  will  also  hold  if  the  curvature  equals 
^  (see  Note  2.2),  provided  certain  other  conditions  hold.  The  reason  is  that  if  the 
curvature  equals  £  the  double-counted  points  have  zero  measure,  since  /j  is  a  measure 
in  3  dimensions.  The  key  is  preventing  Wr(p,  At)  from  doubling  back  on  itself  with 
any  positive  p-?,  measure  (i.e.,  volume).  In  Assumption  2.4,  we  required  that  there 
are  no  intersections  (even  ones  of  zero  measure)  since  the  equality  specified  there 
was  not  in  measure,  but  in  the  actual  spatial  points  themselves. 

We  conclude  this  section  by  explaining  why  the  leading  edge  is  used  as  opposed 
to  the  more  popular  disk  version.  The  reason  lies  in  Lemma  2.3.  Since  we  use  a 
general  density  function  /,  it  is  imperative  that  as  A  moves  forward,  we  add  the 
proper  points,  s  6  S',  to  the  encounter  window.  The  points  to  add  are  those  on  the 
leading  edge.  Most  of  the  literature  in  Kinetic  Theory,  on  the  other  hand,  assumes 
/  is  uniform  in  S,  so  from  Lemma  2.1  it  makes  no  difference  which  points  are  added, 
provided  the  same  volume  is  added.  As  seen  in  Theorem  2.1  the  volume  added  to 
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the  encounter  window  as  A  moves  forward  a  certain  distance  is  equivalent  to 

the  volume  of  a  simple  cylinder.  Therefore  kinetic  theory  prefers  the  more  easily 
visualized  cylinder  formed  by  a  disk  perpendicular  to  the  direction  of  motion  and 
integrated  over  the  path. 

2.3  v  as  a  function  of  S 

The  previous  section  laid  ground  work  by  defining  bodies,  motions,  and  en¬ 
counter  windows.  The  focus  of  this  section  is  to  bring  all  of  these  together  into  a 
single  encounter  model.  Kinetic  theory  provides  a  starting  point  for  our  develop¬ 
ment,  since  much  of  the  work  regarding  relative  velocities  has  already  been  done.  It 
is  generally  known  in  the  held  of  kinetic  theory  that  the  expected  relative  velocity 
between  two  bodies  moving  according  to  the  Assumptions  2.2  and  2.3,  integrated 
over  all  possible  collision  angles  in  3-space  is  given  by  Equation  (2.1) (for  a  proof  see 
Loeb  (1934)  or  Chapter  3  of  this  thesis  ).  This  result  was  first  published  by  Clausius 
(1858)  and  has  been  extended  in  many  ways.  To  the  author’s  knowledge,  however, 
it  has  never  been  extended  as  it  is  in  the  following  lemma  where  we  allow  the  speeds 
themselves  to  be  functionally  related  to  the  space. 

Lemma  2.4  If  A  is  moving  along  pit )  in  S  at  speed  va,  and  B  e  B  is  located  in 
S  and  moves  in  random  directions  with  locally  defined  speeds  vjs{s  €  S)  then  the 
average  relative  speed  between  A  and  B  is 

E  s[Vrel\  =  Eg[vrei(vA,VB(s))}j^^ds  (2.5) 

where  Eg[vrei(vA,  vb)\  is  given  by  Equation  (2.1) 

Proof.  Noting  that  the  probability  that  B  is  located  at  s  is  f(s),  the 
remainder  of  the  proof  follows  from  the  definition  of  the  mean.  ■ 
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From  the  formula  of  Equation  (2.1)  and  Lemma  2.4,  it  can  be  seen  that  the 
results  are  the  same  if  va  varies  over  S  and  vb  is  constant.  The  reason  is  because 
the  speed  of  A  relative  to  B  is  equal  to  the  speed  of  B  relative  to  A  (the  velocity 
simply  points  in  the  opposite  direction).  Thus,  our  assumptions  that  va  is  constant 
and  us  is  a  function  of  s  is  without  loss  of  generality. 

This  shift  in  the  speed  of  the  encountering  particle,  A,  is  of  great  importance, 
because  velocity  is  intimately  related  to  distance  as  seen  in  Definition  2.1.  Since, 
however,  the  actual  distance  covered  by  A  during  some  time  St  is  independent  of 
the  motion  of  B  (only  the  relative  distance),  we  adopt  a  new  scaled  measure  that 
dilates  the  standard  Lebesgue  measure,  /xi,  and  so  accounts  for  the  relative  increase 
in  speed. 

Definition  2.6  Suppose  D  is  a  d-dimensional  cube.  Define  as  ^e  Lebesgue 

measure  of  f!  where  the  first  dimension  is  dilated  by  the  ratio  of  the  relative  speed 
between  va  and  vb  divided  by  va ■  Mathematically, 

(va,vb)  Vrei(v  Aivb)  /n  \  /n  \ 

/fid, re/)  ~  - -=* - /Xi(D{i})/Xd_i(D{2,3...,d}) 

=  Vrel^B\(ti)  (2.6) 

va 

where  Qn  refers  to  the  \n\- dimensional  cube  created  by  the  intersection  of  Q  and  the 
sub-space  of  dimensions  given  in  n. 

The  reason  that  the  dilation  is  only  in  the  first  dimension  is  because  from  Definition 
2.1  v  maps  to  M,  so  the  dilation  due  to  relative  velocity  also  only  occurs  in  M.  LIsing 
the  assumption  that  our  measures  are  in  M3,  we  will  use  the  shorthand  = 


Note  2.3  Suppose  A  and  B  are  moving  in  S  with  speeds  i>a  and  i>b  respectively. 
Then  it  is  easy  to  see  from  Equation  (2.6)  that  though  the  static  encounter  windows  of 
A  and  B  are  different  (assuming  va  vb),  their  encounter  windows  relative  to  each 
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other  are  of  equal  measure.  This  result  leads  to  the  following  intuitive  expectation 
for  relative  encounter  windows. 


Theorem  2.2  The  expected  measure  of  the  new  encounter  window  swept  by  A  rel¬ 
ative  to  a  moving  B  during  an  infinitesimal  time  interval  St  =  (t,  t  +  dt]  C  At, 
conditional  on  the  speeds  of  A  and  B,  is  given  by 

E  =  E^i^a,vb)}  ^w*(p,St))dt 

va 


Proof.  The  proof  follows  from  Equation  (2.6)  and  (2.4) 


E[ddA,l(W/;(P,«))]  =  E 


Vrel(vA,VB) 


va 


p(Wf(p,  St) 


=  E 


Vrel{v Ai  V b)  2 


VA 

E  [Vrel(,VA,VB)\ 
VA 


nr2VA\St\ 


T(Wf(p,  St)) 


Then  by  combining  Lemma  2.1  with  Equation  (2.6),  we  can  show  a  similar 
property  on  the  /  measure  of  Wf(p,  St): 

Mva,vb)(w;(p,  st)  =  tw‘iVA'VB>  f(w;(p,st))  (2.7) 

VA 

Which  presents  us  with  a  dilemma:  frei  is  no  longer  a  probability  measure.  Depend¬ 
ing  on  va  and  vB  and  f(Wf(p,5t)),  it  is  feasible  to  generate  frei  measures  that  are 
larger  than  unity.  Clearly  this  makes  the  use  of  the  binomial  distribution  for  the 
count  process  of  X  infeasible. 
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2.4  The  Extended  Spatial  Poisson  Process 


The  reason  that  frei  can  grow  beyond  unity  is  because  if  motion  is  possible,  it 
is  feasible  for  a  single  element  B  e  B  to  exist  in  two  distinct  locations  in  W*  (p,  St) 
during  the  time  St.  This  would  occur  when  B  was  discovered  at  point  s  at  time  t, 
and  then  it  moved  to  point  s'  where  it  was  discovered  at  time  t' .  However,  if  we 
could  somehow  remove  a  point  from  the  candidate  list  of  encounter-  able  elements 
of  B  immediately  after  it  was  detected,  the  problem  would  be  eliminated.  The 
Poisson  distribution  describes  the  limiting  behavior  required  here,  except  that  it 
never  changes  the  arrival  rate.  Here  we  use  the  extended  Poisson  process  with  a 
linearly  decreasing  arrival  rate  at  each  transition,  first  described  by  Faddy  (1990). 

As  M.  Faddy  (1997)  explains,  a  binomial  distribution  gives  rise  to  the  linearly 
decreasing  sequence  of  arrival  rates  for  a  birth-death  process:  Xn  =  —(N  —  n)  ln(l  — 
p)  for  0  <  n  <  N.  No  provision  has  been  made  in  the  literature  to  allow  for 
inhomogeneity  of  An,  so  we  will  use  the  expectation  of  f(W*(p,St ))  to  give  us  a 
constant  arrival  rate  over  the  incremental  encounter  window. 


symbol  [Arrival  ratejA  =  -  ln(l-E  [f(W*(p,  St))])  =  -  In  ^1  -  (2-8) 

Then  the  rate  transition  matrix  of  the  extended  spatial  Poisson  process  corresponding 
to  the  spatial  binomial  process  X  is  given  by: 


Q 


-NX  NX  0 

0  -(N-1)X  (N-1)X 

0  0  -(N-2)X 


0 

0 

0 


0  0  0  •••  0 


(2.9) 


This  pure  birth  processes  can  be  propagated  in  time  using  the  Chapman-Kolmogorov 
Forward  equations.  In  the  spatial  case,  we  integrate  the  equations  over  S  instead  of 
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t.  Substituting  p(W*(p,5t ))  in  for  t  in  (2.4)  of  Faddy  (1997),  and  setting  the  initial 
condition  p a{W*{p,  St))  =  [1, 0,  0, ... ,  0],  we  will  find  the  probability  distribution  for 
the  number  of  encounters  in  the  region  W* (p,  St),  using  the  formula 

Pn  =  p0exp(Q  *  p(W*  (p,  St)))  (2.10) 

This  can  be  shown  to  be  equal  to  the  binomial  distribution  Binomial  (TV,  f(W*  (p,  St))), 
where  the  nth  element  of  pn  corresponds  to  the  binomial  probability  of  n  successes. 

Extending  Equation  (2.10)  to  the  problem  where  the  elements  of  B  are  moving 
we  can  take  two  equivalent  perspectives  to  account  for  the  relativity.  The  first 
perspective  was  the  one  presented  previously  in  Definition  2.6,  namely,  in  order  to 
require  the  same  time  to  cover  the  same  absolute  distance  at  an  increase  relative 
velocity,  the  Lebesgue  measure  in  the  direction  of  relative  motion  must  be  dilated 
in  proportion  to  the  increase  in  velocity.  This  would  be  equivalent  in  the  extended 
Poisson  process  to  dilating  the  time  t  by  some  proportion.  Time  dilation  is  not 
frequently  studied  in  the  area  of  counting  processes.  However,  the  issue  of  relative 
speeds  and  the  resultant  perception  of  time  or  space  dilation  is  extremely  important 
for  dynamic  spatial  processes. 

If  one  prefers  to  avoid  the  concept  of  space  or  time  dilation,  an  alternative 
perspective  is  that  there  is  an  apparent  increase  in  the  number  of  encounter- able 
elements  in  B.  This  perspective  derives  from  the  fact  that  the  increased  number 
of  collisions  can  come  by  either  increasing  the  density  or  inspecting  a  larger  space. 
Note,  however,  that  the  coefficient  is  applied  only  to  the  remaining  encounter-able 
elements,  thus,  the  coefficient  is  placed  outside  the  parenthesis  in  the  elements  of  Q. 

It  can  be  seen  from  Equation  (2.9)  and  (2.10),  that  the  effect  is  the  same 
regardless  of  which  perspective  is  taken.  The  exponential  term  is  increased  by  the 
coefficient  vrei(vA,  vb),  except  that  in  the  perspective  of  increasing  the  remaining  N, 
the  coefficient  is  buried  in  Q. 
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Example  2.1  Reversing  the  roles  of  A  and  B,  we  consider  the  case  where  the  veloc¬ 
ity  of  A  is  changing  deterministically  from  the  influence  of  constant  uniform  acceler¬ 
ation  and  elements  of  B  are  moving  with  constant  speed  in  random  directions.  Recall 
from  Note  2.3  that  this  formulation  is  equivalent  to  the  one  where  va  is  constant  and 
vb  varies  by  spatial  location. 

LetvA  =  Vo ~gt  where  vq  =  (t>o,0,Vo)  (in  Euclidean  coordinates),  g  =  (0,0,  9.8) 
and  t  G  [0,2^;^].  We  assume  A’s  path  has  the  starting  point  p(  0)  =  (0,0,0).  From 
this  initial  condition  using  the  velocities  from  va,  we  get  the  path 

P(t )  =  {(x,y,z)  :x  =  2t0f1  y  =  0,  2  =  v(t  -  ]-gt2} 

vb  1 

where  h  =  {v(f)2/{2 g)  and  a  =  v^Vy/g.  This  describes  ballistic  motion  along  a 
parabolic  arc  from  the  origin  to  the  point  (2a,  0,  0)  with  zenith  ( a,0,h ).  R  can  be 
shown  that  this  path  meets  the  requirements  of  Assumption  2.f  if  h  >  r. 

The  velocity  of  A  as  a  function  of  z  is  given  by 


va(x,  y,  z)  =  (v%,  0,  (t^)2  -  2 zg) 


(2.11) 


for  ( x,y,z )  :  x  G  [0,2a],  y  =  0,2  =  h  (l  —  j .  The  incremental  encounter 
window  swept  by  A  moving  along  p  has  measure  (Reese  and  Sondow,  2009) 


R(W*) 


nr2  |  -Va2  +  4 h2  +  —  In 


i2 

4  h 


'2h  +  Va2  +  4:h2' 


where  we  use  the  abbreviated  form  W*  =  Wf  (p,  St).  Assuming  that  X  is  uniform  in 
S,  the  probability  that  a  stationary  B  is  in  this  window  is 

To  make  this  example  more  concrete,  suppose  we  start  the  model  with  the  inputs 
shown  in  Table  2.1.  From  the  formulas  above,  we  can  tell  that  the  zenith  occurs  at 
altitude  h  ~  11.5m,  the  total  shot  distance  is  2a  ~  2  *  30.6,  and  the  duration  of 
the  flight  is  tn  ~  4.1  seconds.  With  this  information  we  can  be  confident  that  the 
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Table  2.1  Velocity  inputs  for  example  problem 


Vo 

(15  m/s,  0,  20  m/s) 

9 

9.8  m/s2 

vb 

20  m/s 

Table  2.2  Scenario  inputs  for  example  problem 


S 

[0, 100]  x  [0,  20]  x  [0, 15]  =*►  fi(S)  =  3000m3 

r 

lm 

N 

4 

parameters  shown  in  Table  2.2  for  the  scenario  are  reasonably  spacious  to  encompass 
W* .  They  provide  a  sufficient  space  to  fully  contain  the  collision  window. 

This  gives  us  a  collision  window  volume  of  p(W*)  ~  105m3.  If  we  assume  that 
the  4  elements  of  B  are  uniformly  distributed  in  S  and  stationary,  then  the  density 
function  is  for  a  unit  of  volume 

!  _  ! 

JfS)  ~  3000 

and  from  Equation  (2.8)  we  compute  A  =  —  ln(l  —  f(W*))  ~  0.0012,  since  dF(W*)  = 
7rr2/u(dp)/3000  f(W*)  =  7rr2/3000.  The  resultant  extended  Poisson  process  Q 

matrix  then  becomes: 


-.0042  .0042  0  0  0 

0  -.0031  .0031  0  0 

Q  ~  0  0  -.0021  .0021  0 

0  0  0  -.0010  .0010 

0  0  0  0  0 

Now  if  the  elements  of  B  move  at  speeds  vb  in  random  directions  in  M3,  we  can 
integrate  Equation  (2.5)  to  find  the  mean  relative  velocity  between  A  and  an  arbitrary 
element  B  e  B.  Since  the  integration  is  rather  lengthy,  we  resort  to  numerical 
techniques  to  generate  the  expected  relative  velocity. 
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Numerically  integrating  with  3  point  Gaussian  Quadrature  (Burden  and  Faires, 
2005)  over  all  altitudes  we  find  that  the  expected  relative  velocity  for  this  scenario 
is  approximately  29.92 m/s.  This  means  that  the  average  relative  velocity  between  A 
and  some  arbitrary  element  B  e  B  is  approximately  20%  higher  than  it  would  be  if 
the  B  element  were  stationary.  To  incorporate  the  effect  of  increasing  the  speed  on 
the  collision  probability,  we  can  either  increase  each  (. N  —  i)  in  Q  by  20%  or  we  can 
increase  t  by  the  same  proportion.  Here  we  choose  to  increase  (N  —  i).  This  gives 
us  a  new  Q  matrix: 

-.0050  .0050  0  0  0 

0  -.0038  .0038  0  0 

Q  ~  0  0  -.0025  .0025  0 

0  0  0  -.0013  .0013 

0  0  0  0  0 

Finally,  we  can  exponentiate  Q fi{W*)  to  find  the  transition  probability  matrix 

P 

v 

0.5326  0.3622  0.0940  0.0107  0.0005 

0  0.6195  0.3223  0.0549  0.0033 

P  ij  =  0  0  0.7298  0.2482  0.0220 

0  0  0  0.8489  0.1511 

0  0  0  0  1 

We  use  the  initial  condition  of  0  encounters  with  xq  =  [1,0,  0,0,0]  (alterna¬ 
tively  we  could  use  the  encounter  probability  for  the  ball  around  the  starting  position, 
/3r(0,0,0)),  to  generate  the  probability  distribution  for  the  number  of  encounters  at 
tn;  %t„- 

xtn  =  xoPij  =  [0.5326,  0.3622,  0.0940,  0.0107,  0.0005] 
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2.5  Conclusion 


In  summary,  it  is  clear  that  the  extended  Poisson  process  is  applicable  to  spatial 
point  processes,  and  can  provide  significant  advances  both  from  a  statistical  perspec¬ 
tive,  by  accounting  for  over-  and  under-dispersed  data,  as  well  as  the  probabilistic 
perspective,  in  that  it  unifies  a  variety  of  spatial  point  processes  into  one  common 
framework.  This  framework  also  has  a  unique  ability  to  address  motions  concerns 
at  the  infinitesimal  level,  similar  to  the  spatial  Poisson  process,  while  maintaining 
the  features  of  alternative  counting  processes.  Thus,  in  the  field  of  dynamic  spatial 
processes,  it  opens  up  a  whole  new  horizon  of  feasible  encounter  modeling.  The 
applications  for  such  studies  include:  hard-body  collisions  in  finite  particle  systems, 
search  and  rescue  operations,  and  birth-death  processes  for  links  on  a  mobile  ad-hoc 
wireless  networks. 

On  the  other  hand,  the  methods  currently  in  use  in  spatial  process  models  may 
provide  some  insight  to  advance  the  theory  of  extended  Poisson  processes.  To  date, 
the  processes  are  limited  in  that  the  rates  only  depend  on  the  number  of  arrivals 
at  some  time,  not  the  current  time  itself.  Cox  and  Gibbs  processes  (Mpllcr  and 
Waagepeterson,  2004)  suggest  that  there  may  be  parallel  constructs  in  the  temporal 
domain  to  allow  for  time- variant  versions  of  the  extended  Poisson  process  model  with 
time-intransient  Markov  chains. 
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3.  Case  Study 

Assessing  Probability  of  Aircraft  Fratricide 
Due  to  Artillery  Fires 


3.1  Introduction  and  Problem  Description 

The  risk  analysis  for  artillery  fires  colliding  with  aircraft  in  a  conflicted  airspace 
is  increasingly  becoming  a  source  of  concern  for  mission  planners,  decision  makers, 
and  commanders.  In  December  2007,  the  Center  for  Army  Lessons  Learned  pub¬ 
lished  a  requirement  for  an  analytical  model  for  assessing  the  risk  of  aircraft  fratricide 
due  to  fires  (Neuenswander,  2008).  Research  conducted  by  the  Air  Force  Doctrine 
Development  and  Education  Center  (Bethea  and  Herbranson)  introduces  a  probabil¬ 
ity  model  that  introduces  the  concept  of  a  shell  creating  a  collision  “tube”  through 
the  sky  with  a  geometric  approximation  of  the  ballistic  path.  Their  model,  however, 
fails  to  account  for  the  dynamic  motions  of  the  aircraft  and  artillery  projectiles. 
The  National  Defense  Industrial  Association  (Strike,  Land  Attack  and  Air  Defense 
(SLAAD)  Committee  of  the  National  Defense  Industrial  Association  (NDIA),  2001) 
also  conducts  a  simulation  study  that  concluded  a  1%  risk  of  fratricide.  Despite  the 
increase  in  studies  on  the  subject,  no  published  model  exists  to  date  for  probabilistic 
estimation  of  the  risk  of  fratricide  of  fixed  wing  aircraft  due  to  artillery  fires. 

Figure  3.1  shows  a  picture  describing  the  problem  scenario.  Artillery  units 
fire  indirectly  (ie.  the  gunner  cannot  see  the  target)  arcing  through  an  altitude 
range  where  aircraft  are  flying.  The  model  must  take  into  account  all  the  important 
parameters,  but  be  flexible  so  that  it  can  operate  with  minimal  information. 

The  approach  to  solving  this  problem  begins  in  Section  3.2  with  a  literature 
review  on  general  collision  models.  The  best  candidate  for  solving  this  problem  is 
selected  and  adjusted  to  better  match  the  particular  scenario  given  by  this  problem. 
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Aircraft 

altitude 


Figure  1  -  Problem  scenario 


Figure  3.1  The  Problem  Scenario 

Specifically,  the  model  is  adjusted  in  Section  3.3.1  to  account  for  ballistic  trajectory, 
spatial  M/D/oo  birth-death  processes  in  Section  3.3.2,  and  Section  3.3.3  accounts  for 
uncertainty  in  the  parameters.  After  these  adjustments,  the  final  model  represents  an 
original  contribution  to  the  literature  of  collision  models.  The  model  is  implemented 
in  software  and  subjected  to  sensitivity  analyses.  Finally,  numerical  results  are 
computed  for  a  simple  scenario  and  some  conclusions  are  drawn. 


3.2  Finding  the  Right  Model 

The  issue  of  things  colliding  is  a  point  of  interest  in  many  different  fields  of 
study.  The  literature  is  replete  with  collision  studies  that  represent  many  different 
physical  collision  models  with  different  mathematical  bases.  Most  of  these  studies 
seem  to  solve  their  respective  problem  in  one  of  three  methods  briefly  presented  here. 
For  further  discussion  on  the  models  not  carried  through  in  this  study,  see  Chapter 
1. 


1.  The  Historical  Model :  The  general  concept  is  to  take  historical  collision  data 
from  similar  scenarios,  and  use  these  to  estimate  future  collision  probabilities. 
It  is  a  very  flexible  method  that  is  simultaneously  very  simple.  The  draw¬ 
backs  are  finding  the  data  or  developing  and  validating  a  realistic  simulation 
to  generate  the  data. 
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2.  The  Gas  Particle  Model.  Here  the  incident  particle  moves  through  space, 
sweeping  through  a  collision  window  in  the  process.  Spatial  point  processes 
(SPP)  provide  a  probability  of  a  target  particle  existing  in  the  window.  The 
strength  of  the  model  is  its  wide  applicability.  Its  primary  weakness  is  that 
it  quickly  becomes  unwieldy  and  complex  when  assumptions  of  independence 
and  uniformity  are  lost. 

3.  The  Satellite  Model  This  third  model  uses  geometric  techniques  on  location 
and  path  information  of  all  objects  to  compute  the  points  of  closest  approach 
for  every  incident-target  object  pair.  Then  it  uses  a  Gaussian  error  distribution 
at  the  point  of  closest  approach  to  estimate  an  upper  bound  on  the  probability 
of  collision.  The  strength  of  this  model  is  its  precision,  but  its  weakness  is  the 
large  data  requirements  and  the  probability  assumption  of  normality. 

As  we  consider  these  three  models  for  use  in  the  problem  at  hand,  we  must  eval¬ 
uate  their  applicability  to  the  problem,  flexibility  to  work  under  uncertainty,  and 
simplicity  of  implementation.  We  quickly  reject  the  satellite  model  for  its  data  re¬ 
quirements,  which  greatly  exceeds  the  amount  of  data  that  is  generally  available  in 
real-time  combat  situations.  This  presents  us  with  a  choice  between  the  historical 
model  and  the  gas  particle  model. 

The  gas  particle  and  the  historical  models  have  both  been  used  to  develop 
midair  collision  estimates.  Alexander  (1970)  uses  the  gas  particle  model  in  estimating 
the  risk  of  mid-air  collisions  for  aircraft  using  the  “free  flight”  rules.  His  “back  of 
the  envelope”  calculations  demonstrate  that  the  gas  particle  is  a  viable  method 
of  estimating  random  collisions  between  particles  moving  in  uncertain  directions. 
The  historical  model  was  implemented  using  simulations  by  the  National  Defense 
Industrial  Association  (Strike,  Land  Attack  and  Air  Defense  (SLAAD)  Committee 
of  the  National  Defense  Industrial  Association  (NDIA),  2001)  to  test  the  “Big  Sky 
-  Little  Missile”  theory  for  Joint  Fire  doctrine  development.  They  do  not  provide 
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details  on  the  half-width  confidence  interval  results,  but  report  that  approximately 
1%  of  their  simulated  scenarios  resulted  in  a  fratricide. 

The  historical  model  relies  on  many  past  observations  of  event  phenomenon. 
This  technique  is  not  wise  to  implement  in  combat  theater  models  apart  from  simu¬ 
lation.  However,  in  view  of  the  fact  that  the  estimated  collision  frequency  is  low,  it  is 
reasonable  that  a  simulation  approach  would  require  an  extensively  large  simulation 
time  or  number  of  replications  to  develop  sufficiently  bounded  confidence  intervals 
for  realistic  results.  Nevertheless,  Alexander  demonstrated  that  the  gas  model  can 
be  applied  to  the  problem  at  hand  without  undue  difficulty. 

3.3  Model  Development 

The  model  is  most  easily  understood  in  the  simpler  two  dimensional  case  pre¬ 
sented  by  Kim  (2002).  Kim  examines  the  probability  of  activating  a  mine  while 
traversing  a  minefield.  If  a  person  travels  a  distance  d  across  the  minefield,  and  he 
will  activate  any  mines  within  a  radius  r,  any  mines  located  inside  an  area  of  dr  will 
be  activated.  If  the  mine  density  is  a  mines  per  unit  area,  then  the  expected  num¬ 
ber  of  mines  activated  will  be  dro.  The  aircraft  collision  model  is  in  3  dimensions 
and  the  artillery  shells  are  in  motion,  so  modifications  must  be  made,  but  the  key 
concepts  are  the  same.  Accounting  for  the  third  dimension,  the  expected  number  of 
static  fratricides  in  a  collision  window  swept  by  the  aircraft  is 

nr2  do  (3.1) 

Visually  the  model  we  use  here  is  presented  in  Figure  3.2.  In  Section  3.3.1  we 
use  relative  velocities  to  add  in  motion  of  the  shells.  The  second  step,  described  in 
Section  3.3.2,  is  to  estimate  the  shell  density  in  the  collision  window  by  means  of 
a  spatial  point  counting  process.  Finally,  these  two  results  are  combined  in  Section 

3.3.3  to  give  a  numerical  estimate  of  the  collision  risk  while  accounting  for  uncertainty 
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Figure  3.2  Model  flow  chart. 


in  the  parameters.  In  forms  of  calculations  similar  to  the  formulas  of  Kim’s  model, 
the  distance  d  is  found  in  Step  1  of  the  model  and  a  is  calculated  in  Step  2. 

3.3.1  Step  1:  Relative  Velocity  (Applying  Ballistic  Formulas) 

This  collision  model  deviates  from  Kim’s  minefield  model  in  that  the  shells  are 
in  motion  whereas  mines  are  not.  Fortunately,  Loeb  (1934)  provides  a  method  to 
account  for  shell  motion  using  relative  speeds  between  the  aircraft,  va,  and  the  shell, 
vs.  Mathematically,  this  is  computed  using  the  Law  of  Cosines  (see  Figure  3.3): 

Vrei  =  c  =  V a2  +  b2  —  2 ab  cos  9  =  \J v\  +  v2s  —  2vavs  cos  d 

In  our  problem  a  and  b  are  the  speeds  va  and  vs.  d  is  the  angle  of  incidence 
between  the  two  motion  vectors.  The  expected  relative  speed  is  found  by  integrating 
over  all  possible  d  angles  as  follows.  We  assume  that  the  aircraft  path  is  in  some  given 
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Figure  3.3  The  geometry  of  relative  velocities. 


direction  and  that  the  artillery  shell  could  be  coming  at  it  from  any  direction.  Thus, 
without  loss  of  generality,  we  fix  the  aircraft’s  direction  of  motion  to  the  vertical 
axis:  va  =  z.  Then  we  assume  that  the  shell’s  motion  vector  is  uniform  in  3-space. 
If  we  desire  to  test  whether  the  angle  of  approach,  Q  is  between  some  angle  d  and 
-d  +  d-d,  we  must  take  the  ratio  of  the  surface  area  of  the  “belt”  depicted  in  Figure 
3.4  to  the  surface  area  of  the  sphere.  If  we  consider  the  limiting  case  as  d-d  —>  0, 


Figure  3.4  Computing  the  probability  of  angle  of  incidence. 


the  “belt”  approaches  the  side  of  a  cylinder  with  radius  sin  d  and  height  d-d.  So  the 
ratio  is  computed  as  follows: 


/(-d)  =  P  ('&<(<$  +  d-d) 


2n  sin  -dd-d 
47T 


-  sin  -dd-d 
2 


(3.2) 
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We  then  use  the  probability  of  the  incidence  angle  t)  to  find  the  expected 


relative  speed  between  two  bodies  moving  in  3  space. 


E  (vrei)  =  /  Vrei(va,  Vs,  <(<'#  +  dl?) 


'0=0 


'0= 0 

1 


yj v\  +  f  2  —  2vavs  cos 'd  (  -  sin  i)<\i) 


6  VnV 


{y2s  +  v2  —  2vav2s  cos '()) 


2/3' 


a^s  u 

l 


J  0=0 


W  +  vs  ~  2vavs)3/2  +  (v2a  +  v2s  +  2nans) 


3/2' 


cms  ■- 

1 


(( 


,2\3/2 


aus  *- 

1  r 


6vnvq 


Va  -fS||J  )  +((W+^s)) 

3  +  (fa  +  fs)3] 


2\3/2 


fa  -  V. 


Vs+  3/y,  vs  >  va 


va+  fS<fa 


(3.3) 


Note  that  in  Equation  (3.3),  the  speeds  are  considered  constant.  However,  due 
to  the  physics  of  ballistic  motion,  shell  velocity,  is  not  constant  over  the  shell’s 
arc.  Rather  the  shell  speed  decreases  as  it  goes  up  and  increases  again  as  it  comes 
down.  In  general,  vs  is  a  function  of  launch  angle,  9;  muzzle  velocity,  vm ;  vertical 
position,  y;  and  acceleration  clue  to  gravity,  the  constant  g.  As  an  initial  model  we 
use  the  ideal  ballistic  arc,  wherein  the  instantaneous  shell  speed  at  any  given  altitude 
follows  the  equation  (for  a  proof,  see  Chapter  2) 


v,{6,y,vm)  =  ■/[«„, cos#]2  +  [«,„siii0]'2  -  2 yg 


(3.4) 


where  9  is  the  angle  of  elevation  at  launch,  vm  is  the  muzzle  velocity  at  launch,  g  is 
the  gravitational  constant  (32.2  ms-2),  and  y  is  altitude. 

Then  we  set  W  to  be  the  window  volume  swept  by  the  aircraft  during  the 
conflicted  time.  If  we  didn’t  have  relative  velocity  to  consider,  we  would  have  W  = 
7r r2tva.  However,  we  do  take  the  relative  velocity  into  account  using  Equation  (3.3). 
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Moreover,  we  also  consider  the  speed  of  the  aircraft  at  a  particular  height.  Given 
the  launch  angle,  6  and  the  muzzle  velocity  vm,  we  combine  Equations  (3.1),  (3.3), 
and  (3.4)  to  find  the  new  collision  window  formula: 

pin 

Wrd=  irr2E{vrei(va,vs(y(t)),tf)}dt  (3.5) 

Jto 

where  [t0,tn\  is  the  time  that  the  shell  spends  in  the  conflicted  airspace. 

3.3.2  Step  2:  Shell  Density  (Applying  Spatial  Point  Process  Theory) 

Given  the  volume  Wrei,  we  next  count  how  many  shells  are  actually  in  the 
window.  The  appropriate  tool  for  such  a  counting  process  is  called  a  spatial  point 
process  (SPP).  SPP  theory  is  primarily  used  for  statistical  analysis  used  in  areas 
such  as  forestry,  ecology,  seismology,  epidemiology  and  material  science  (Mpller  and 
Waagepeterson,  2006).  SPP’s  come  in  different  variations  based  on  how  the  points 
are  distributed  in  space.  The  Binomial  and  Poisson  distributions  are  frequently 
used  to  model  data  that  is  fairly  evenly  spread  out  throughout  the  space.  More 
advanced  models  such  as  Cox  and  Gibbs  processes  are  used  if  the  points  tend  to 
cluster  together  or  scatter  apart.  Rarely  do  the  objects  move  in  any  of  the  models 
except  in  two  domains:  kinetic  theory  (Brush,  1965)  and  orbital  debris  environments 
(Jenkin  and  Gick,  2002).  In  both  of  these  cases,  the  Poisson  SPP  is  the  model  of 
choice  because  of  its  memoryless  property  which  creates  truly  random  distributions. 
In  the  space/time  context  the  memoryless  property  is  as  follows:  as  I  fly  through 
the  space,  assuming  that  I  do  not  double  back  on  myself,  the  number  of  collisions 
I’ve  had  so  far  has  no  bearing  on  the  number  of  collisions  I  will  have  in  the  future. 

We  assume  that  the  Poisson  process  is  a  reasonable  estimate  for  the  shells 
entering  the  airspace.  This  corresponds  to  an  exponentially  distributed  inter-firing 
time  for  the  artillery  batteries.  The  assumption  of  a  exponential  distribution  may 
seem  rather  arbitrary,  but  it  is  justified  on  the  fact  that  the  distribution  is  the 
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most  random  law  that  can  be  used,  and  so  provides  a  conservative  bound  on  the 
risk.  Furthermore,  as  the  number  of  artillery  units  firing  into  the  space  increases, 
Khintchine’s  (Gross  and  Harris,  1998)  law  dictates  that  the  inter-firing  times  will 
converge  to  an  exponential  distribution. 

However,  we  must  also  take  into  account  the  fact  that  the  shells  fly  back  out  of 
the  airspace.  This  requires  using  a  spatial  birth-death  point  process.  Since  we  know 
from  the  ballistic  arc  exactly  how  long  the  shell  remains  in  the  altitude  window,  we 
use  a  deterministic  formula.  If  an  aircraft  is  flying  between  altitudes  ymm  and  ymax , 
then  the  duration  that  the  shell  spends  inside  the  range  [ymm,  ymax]  is 


2 g  min (amin, 


h)  - 


2 g  min (amax, 


(3.6) 


where  vv  is  the  vertical  component  of  the  muzzle  velocity  (vv  =  vmsin#),  h=  v^/2 g  is 
the  zenith  (henceforth,  max  ord)  of  the  ballistic  arc  (again  assuming  an  ideal  ballistic 
arc),  and  the  construct  armn/rnax  =  mm{ymin/max,  h}  handles  the  cases  where  the  shell 
falls  short  of  the  min  or  max  altitude  of  the  aircraft. 

In  summary  we  assume  that  shells  are  fired  into  the  space  at  exponentially 
distributed  time  intervals  with  an  average  rate  A  and  that  each  shell  spends  the 
same  amount  of  time  in  the  collision  window  given  by  t  in  Equation  (3.6).  In  queuing 
theory  notation,  this  is  represented  as  an  M/D/oo  birth-death  process.  Adjusting 
the  formulas  of  Gross  and  Harris  from  the  M/D/c  queue  (Gross  and  Harris,  1998), 
the  number  of  shells  in  the  conflicted  airspace  volume  at  any  point  in  time  is 
Poisson[Af]  (~  means  “distributed  as”).  Thus  the  number  of  shells  in  the  whole 
collision  altitude  at  any  given  time  is  found  by  counting  how  many  shells  were  fired 
in  the  last  t  seconds.  Since  L  is  the  number  of  shells  in  the  whole  space,  assuming 
that  the  shells  are  uniformly  distributed  throughout  the  space,  we  find  the  shell 
density  for  a  unit  of  air  volume,  a  by  dividing  L  by  the  total  volume  of  conflicted 
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airspace  Wt'- 


L 


rv-/  Poisson 


A  t 

Wr 


(3.7) 


Another  nice  feature  of  the  birth-death  process  as  opposed  to  other  SPP  models 
is  that  it  changes  the  requirements  for  the  memoryless  property.  As  stated  above,  the 
memoryless  property  requires  that  the  aircraft  cannot  double  back  on  itself.  Clearly 
this  presents  a  difficulty  for  estimating  fratricide  risks  for  orbiting  aircraft.  However, 
with  the  spatial  birth  death  process,  we  only  require  that  the  aircraft  cannot  go 
over  its  own  tracks  inside  oft  seconds.  Depending  on  the  altitude  window  that  the 
aircraft  is  flying  in,  t  is  usually  no  more  than  only  a  few  seconds  -  well  inside  of  the 
time  it  takes  for  any  aircraft  to  complete  one  turn  of  an  orbit. 


3.3.3  Step  3:  Probabilistic  Estimation  (Applying  Numerical  Integration) 

We  are  now  able  to  give  a  probabilistic  estimate  of  the  number  of  fratricides 
using  Equation  (3.1).  We  substitute  Equation  (3.5)  in  for  irr2d,  and  Equation  (3.7) 
in  for  a  to  get: 


A  frats 


Nx  =  <7  x  *  Wrei  rs-/  Poisson 


Wrel 

WT 


A  t 


(3.8) 


However,  looking  back  at  Equation  (3.5),  we  see  that  Wrei  depends  on  the 
altitude,  y.  Similarly  t  from  Equation  (3.6)  depends  on  the  muzzle  velocity,  vm ;  the 
angle  of  fire,  9;  and  the  aircraft  min/max  altitudes,  ymtn  and  ymax-  In  addition,  a 
depends  on  Wt  in  Equation  (3.7).  Rarely,  however,  will  the  input  parameters  be 
known  with  certainty.  It  is  for  that  reason  that  Va(jj)  is  left  in  its  functional  form, 
since  rarely  will  aircraft  fly  at  one  specific  altitude.  Similarly,  it  may  be  reasonable 
to  expect  uncertainty  in  other  parameters  as  well.  We  expect  that  Wt,  ymin  and 
Umax  WH1  n°t  change  over  a  given  altitude,  but  it  is  certainly  reasonable  to  expect 
some  variance  in  vm,  9 ,  and  y. 
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Considering  this,  we  parameterize  the  Poisson  distribution  in  Equation  (3.8) 


with 


A*  =  ^Ai 
XV t 


(3,9) 


and  then  we  make  A*  a  function  of  vm,  9 ,  and  y.  Thus  Equation  (3.8)  becomes 


#  frats  =  Nx  =  ax*  Wrei  ~  Poi  [A  *  (y,  9 ,  vm) 


(3.10) 


Then  we  allow  any  probability  distributions  fVm,  fg  and  fy  for  their  respective 
variables.  For  example,  we  will  here  use  uniform  distributions  on  the  launch  angle, 
6  ~  Uni (9min,9max);  muzzle  velocity,  vrn  ~  \Jm([vm]min,  [um]maa;);  and  collision 
altitude,  y  ~  Vm(ymin,ymax)  This  allows  us  to  compute  the  probability  for  a  certain 
number  of  fratricides  by  a  conditioning  argument 


P (Nx  =  n)  =  /  P (Nx  =  n\y,  9,  vm)P(y,  6 ,  vm)  d (y,  9,  vm) 

J  (y,0,Vm) 

( ' y max  i(ym)max) 


P (Nx  =  n\y,  9 ,  vm)f(y)f(6)f(vm)dy  d 9  dvm  (3.11) 


{yminfl  liym)  min) 


where  we  know  from  the  Poisson  distribution  of  Nx  that 


P(NX  =  n\y,  9,vm)  = 


e-r(.y,o,vm)  (A *(y,9,vm)y 


n\ 


Since  the  integral  is  more  complex  in  this  case  (considering  the  amount  of 
computation  necessary  to  find  A*),  we  solve  it  numerically  with  3-point  Gaussian 
quadrature.  Higher  order  quadratures  provide  better  estimates,  but  since  multi¬ 
ple  integrals  are  required  the  number  of  function  evaluations  grows  exponentially. 
Therefore  the  lower  order  quadrature  is  chosen  for  simplicity  in  this  proof-of-concept. 

The  quadrature  technique  requires  us  to  build  a  transformation  function  t  from 
[1,  —1]  to  the  variable  domains.  So,  for  a  general  variable  x ,  we  construct  the  function 
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<t>x{tx)  =  ts\(xmax  -  Xmin )  +  \{xmax  +  xmin)  where  tx  e  [-1, 1].  Then  (Burden  and 
Fa.ires,  2005)  demonstrates  that  in  the  general  case, 

ff(x)ix  =  fj(Ut,fmax  ~  Xmmd  tx 

We  note  that  the  term  xmax  —  xmin  will  be  canceled  out  by  the  fx  terms  of  Equation 
(3.11).  After  these  variable  shifts,  (3.11)  becomes: 


(ty,tg,tVrn)=(  1,1,1) 


P(ffx  =  n)  =  - 


P(iV  =  n\<j>y(ty),  (j>g(tg)l<j>Vrn(tVm))dtVlndtodty 


,-l  ,-l) 


(3.12) 


which  is  solved  using  Gaussian  quadrature  by: 


P(NX  =  n ) 


1 

8 


E 

(ty,tg,tvm)e  {-1,0.1}3 


Cy(ty)cg(tg)cvm  {pVrri) 

P(A^  =  n\<l>y(ty),<l>o(te),<l>vm(tvm)) 


(3.13) 


where 

(+,_\  5/9,  *«€{- 1,1} 

Cx(tx)  \ 

^  8/9,  4  =  0 

We  recognize  that  Equation  (3.13)  is  daunting  to  express  in  a  single  formula, 
but  it  is  quite  easy  to  set  up  using  standard  spreadsheet  software  such  as  Microsoft 
®  Office  Excel  ®  (Microsoft®  Corporation,  2006).  This  is  the  approach  we  used  to 
perform  the  sensitivity  analysis  and  to  derive  the  numerical  results  of  the  following 
sections. 


3-4  Sensitivity  Analysis 

Sensitivity  analysis  was  performed  on  the  above  model.  Figure  3.5  below  pro¬ 
vides  a  screenshot  of  the  model  implementation  in  Microsoft®  Office  Excel®.  The 
data  was  generally  designed  to  both  represent  realistic  data  gathered  from  Field 
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Manual  6-40  (HQ  ,  1999)  and  Jane’s  military  resources  (Daly,  2008),  and  to  provide 
insight  into  the  functional  relationships  represented  in  the  equations.  The  aircraft 
speeds  and  altitudes  are  generally  representative  of  unmanned  aerial  vehicles  (UAVs) 
and  the  muzzle  velocities  are  reduced  to  highlight  functional  relationships.  Aircraft 
altitude  was  represented  using  a  uniform  distribution  with  a  spread  of  1000  ft.  All 
other  parameters  are  held  constant  for  each  observation,  and  those  constants  are 
varied  over  many  observations  to  generate  the  output  graphs  shown. 

Comparing  muzzle  velocity  to  collision  probabilities,  until  the  muzzle  velocity 
attains  a  critical  value,  the  projectile  falls  short  of  the  minimum  aircraft  altitude, 
and  the  collision  probability  is  zero.  As  the  max  ord  gradually  travels  through  the 
altitude  range  the  collision  probability  increases  until  it  reaches  the  maximal  value 
when  the  max  ord  equals  the  maximum  altitude  for  the  aircraft.  Then  the  collision 
probability  gradually  decreases  and  levels  off  at  some  constant  value.  The  leveling 
convergence  at  higher  muzzle  velocities  represents  the  ballistic  path  through  the  at- 
risk  altitudes  becoming  straighter.  As  the  muzzle  velocity  increases,  the  probability 
converges  to  the  limiting  case  of  the  projectile  flying  through  the  conflicted  altitude 
in  a  straight  line  at  an  angle  equal  to  the  angle  of  fire.  This  progression  can  be 
seen  in  Figure  3.6.  Since  the  numerical  results  match  logical  outcomes,  this  serves 
as  validation  for  the  model. 

The  significance  of  this  relationship  between  max  ord,  aircraft  altitude,  and 
fratricide  risk  cannot  be  overstated.  It  shows  that  the  point  of  highest  risk  is  given 
by  a  simple  linear  relationship  between  where  the  aircraft  is  flying  and  the  max  ord 
of  the  projectile.  If  the  aircraft  altitude  is  x  ±  d,  and  the  projectile  is  fired  with  a 
max  ord  of  h,  then  the  probability  of  a  fratricide  is  highest  for  aircraft  flying  with 
an  average  altitude  of  x  =  h  —  d. 

The  significance  of  relative  velocities  is  secondary  to  the  muzzle  velocity  at¬ 
taining  this  peak  value.  This  is  displayed  visually  in  Figure  3.7.  The  fratricide 
probability  peaks  when  the  velocity  attains  the  appropriate  level  for  the  aircraft 
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Inputs 

Artillery  Unit 

Muzzle  Velocity 

Maximum 

733 

ft/sec  ^ 

Minimum 

732 

ft/sec 

Angle  of  Fire  (0  =  horz,  90  =  vert) 

Maximum 

41 

0  </ 

Minimum 

40 

0 

Rounds  per  minute 

2 

rpm 

Aircraft 

Speed  (ground  speed) 

40 

ft/sec 

Altitude  (Above  Ground  Level) 

Maximum 

101 

ft  </ 

Minimum 

100 

ft 

Number  of  aircraft 

5 

✓ 

Scenario 

Surface  Area 

4 

nm! 

Min  desired  miss  distance 

100 

ft 

Vulnerable  time  window 

10 

min 

Outputs 

P(frat)  (per  aircraft) 

197/313 

(=  62.94%) 

Collision  rate 

0.4963045 

f  rats/m  in 

E(#  frats) 

4.9630448 

frats 

Var(#  frats) 

4.9630448 

Mean  ±  1  St.  Dev  (#  frats  during  time  window) 

[0,  0.877] 

frats 

N=!  0  P(0  frats )= 

1/143 

(=  0.7%) 

Notes 

None 

Figure  3.5  Three  groups  of  input  (artillery,  aircraft,  and  scenario)  are  required  to 
produce  an  estimate  of  fratricide  along  with  an  uncertainty  measure  of 
the  estimate. 
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Figure  3.6  p  by  vm  and  y. 


Probability  of  Collision  by  Muzzle  Velocity  and  Aircraft 
Altitude  (Aircraft  So^ed:  200  ft/sec) 


0.50%  \ 
0.00% 


840Aircraft  Altitude 
(+/-  250  ft) 


Muzzle  Velocity  (ft/sec) 


4  — « 

s‘I/  — 

®  -uu  f  Altitude 


4000 
E  u  3000 
2  1  2000 
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2000  4000 
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altitude  (max  ord  =  d  +  x  =  250  +  1000  =  1250  ft,  so  critical  muzzle  velocity  is 
approximately  408  ft/sec).  Until  both  the  muzzle  velocity  and  aircraft  speed  reach 
this  critical  velocity,  changes  in  aircraft  speed  minimally  affect  collision  probabili¬ 
ties.  As  the  aircraft  speed  exceeds  the  critical  velocity,  it  begins  to  dominate  the 
equations,  so  that  further  increases  in  aircraft  speed  drive  the  probability  of  collision 
even  higher.  Figure  3.8  demonstrates  that  even  the  effect  of  increasing  the  aircraft 
speed  can  be  dwarfed  by  relatively  small  changes  in  the  collision  radius.  For  exam¬ 
ple,  increasing  the  miss  distance  by  only  2  ft  has  about  the  same  impact  as  raising 
the  aircraft  speed  by  50  miles  per  hour. 

3.5  Numerical  Results 

The  response  values  from  the  model  indicate  surprisingly  high  results.  A  ’’Big 
Sky  -  Little  Missile”  concept  implies  miniscule  collision  probabilities.  However,  in 
very  normal  scenarios  here,  it  is  not  uncommon  for  the  model  to  return  fratricide 
risks  on  the  order  of  1%  to  5%.  With  input  from  an  experienced  Air  Liaison  Officer 
(Glen  Shilland,  2009),  we  generated  three  scenarios  presented  below.  All  three  have 
artillery  fire  with  muzzle  velocities  uniformly  distributed  between  1525  and  1725 
ft/sec  and  angles  of  fire  uniformly  distributed  between  35  and  45  degrees.  These 
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Figure  3.7  p  by  vm  and  va. 


Probability  of  Collision  by  Muzzle  Velocity  and  Aircraft 
Speed  (Aircraftaltitude:  1000) 


Figure  3.8  p  by  va  and  r. 

Probability  of  Collision  by  Miss  Distance  and  Aircraft  Speed 
(Aircraftaltitude:  1000,  Muzzle  Velocity:  400) 


g55  Miss  Distance  (ft) 


Aircraft  spd  (ft/sec) 
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ranges  correspond  to  a  minimum  distance  to  target  of  about  20  km  and  a  maximum 
of  about  28  km.  The  targets  are  uniformly  distributed  throughout  surfaces  areas  (ie. 
a  top  down  view)  called  grid  squares.  Each  grid  square  is  1  km2.  Since  the  target 
areas  are  so  much  smaller  than  the  distances  from  the  artillery  to  the  target,  we 
assume  that  the  shells  only  pass  through  the  altitudes  on  descent,  so  t  in  Equation 
(3.6)  does  not  have  the  2  coefficient. 

3.5.1  Scenario  1 

Our  first  scenario  was  intended  to  be  conservative.  The  artillery  is  firing  at  a 
target  in  1  grid  square  at  a  rate  of  2  rounds  per  minute  (rpm).  Meanwhile  a  2-ship 
flight  of  A-lOs  attack  the  target.  The  A-lOs  fly  at  330  knots  (kts)  between  4500 
and  5500  ft  above  ground  level  (AGL),  proceed  directly  above  their  target,  and  fly 
out  of  the  grid  square.  The  entire  duration  of  the  conflicted  airspace  is  6  seconds 
(the  approximate  time  it  takes  to  cross  a  grid  square  at  330  kts).  We  desire  the 
probability  that  any  rounds  pass  within  50  ft  of  the  center  of  either  aircraft. 

The  result  of  the  scenario  is  that  there  is  less  than  0.001%  risk  to  each  aircraft. 
Overall,  the  probability  that  no  shells  pass  within  50  ft  of  either  aircraft  is  approxi¬ 
mately  99.95%,  or  alternatively  stated,  there  is  a  0.05%  risk  that  a  round  will  pass 
within  50  ft  of  at  least  one  aircraft.  This  is  a  very  low  risk,  however,  the  results  can 
be  very  sensitive  to  the  radius  of  collision.  Figure  3.9  displays  the  effect  of  increasing 
the  minimum  desired  miss  distance  on  the  probability  of  a  shell  passing  within  the 
radius  of  either  aircraft. 

3.5.2  Scenario  2 

The  second  scenario  considers  the  same  artillery  barrage  on  the  target  inside 
of  1  grid  square.  This  time,  2  RQ-7A  Shadow  unmanned  aerial  vehicles  (UAVs)  are 
in  orbit  gathering  intelligence  on  the  battle  development.  They  are  flying  at  75  mph 
between  4800  and  5200  ft  AGL.  The  entire  scenario  lasts  10  minutes,  and  it  is  again 
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Figure  3.9  Comparing  the  probability  of  no  fratricide  to  collision  radii 


Figure  3.10  Comparing  the  probability  of  no  fratricide  to  scenario  duration 

desired  to  know  the  probability  that  any  rounds  pass  within  50  ft  of  the  center  of 
either  aircraft. 

This  scenario  results  in  a  much  higher  risk:  a  4.87%  risk  that  at  least  one  shell 
will  pass  within  50  ft  of  one  of  the  UAVs.  Each  UAV  individually  has  a  2.46%  risk. 
The  major  reason  for  the  increase  in  risk  from  the  last  scenario  and  this  one  is  the 
duration  of  the  airspace  confliction.  We  show  the  effect  that  duration  has  on  the 
risk  in  Figure  3.10. 

3. 5. 3  Scenario  3 

Our  last  scenario  was  intended  to  examine  a  more  congested  airspace.  This 
scenario  increases  the  bring  rate  to  5  rpm,  but  also  spreads  it  out  over  4  grid  squares. 
The  result  of  this  is  that  the  density  is  reduced  by  .75  rounds  per  grid  square  per 
minute.  In  this  scenario  2  helicopters  are  orbiting  the  area  between  1900  and  2100 
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ft  AGL  at  100  kts,  2  UAVs  are  on  surveillance  flying  at  55  kts  and  4800  to  5200  ft 
AGL,  and  2  A-lOs  are  on  orbit  for  CAS  support  at  9500  to  10000  ft  AGL  and  250 
kts.  The  minimum  desired  miss  distance  for  all  aircraft  is  again  50  ft. 

Individually,  the  aircraft  have  fairly  decent  probabilities  of  no  fratricides:  the 
helicopters  have  a  99.51%  chance  of  no  shells,  the  UAVs  have  98.40%  chance,  and 
the  A-10  planes  have  a  probability  of  97.94%  chance  of  both  being  at  least  50  ft  from 
any  shells.  However,  collectively,  there  is  a  5.06%  chance  that  the  center  of  at  least 
one  of  these  aircraft  will  be  with  in  50  ft  of  some  shell. 

3.5.4  Further  discussion 

These  probabilities  are  of  sufficient  risk  to  warrant  attention.  In  comparison, 
if  a  0.1%  risk  of  incapacitation  is  presented  to  ground  troops  supported  by  artillery 
fire  or  close  air  support,  the  situation  is  termed  “Danger  Close,”  which  requires  that 
the  “supported  commander  must  accept  responsibility  for  the  risk  to  friendly  forces 
when  targets  are  inside  the  [danger  close  zone].”  (Joi,  2003).  In  the  second  and  third 
scenarios,  the  risks  are  50  times  as  large  as  the  “Danger  Close”  threshold. 

Furthermore,  when  considering  the  expected  costs  of  flying  through  conflicted 
airspaces,  the  results  are  frequently  non-trivial.  Not  considering  the  probability  of 
kill  (j9fc)  values,  the  0.05%  risk  to  an  A-10,  priced  at  12  million,  corresponds  to  an 
expected  cost  of  the  mission  of  $60,000.  These  are  just  the  “nuts  and  bolts”  costs 
and  do  not  include  effects  such  as  unit  morale  degradation,  loss  of  public  opinion, 
reduced  of  mission  effectiveness  until  repairs,  or  medical  costs.  Most  of  such  missions 
through  conflicted  airspaces  will  not  realize  such  costs,  since  the  probability  of  safe 
passages  are  so  high,  but  it  would  not  be  wise  to  ignore  such  non-trivial  risks  with 
such  large  consequences. 

These  numerical  results  suggest  that  the  current  military  policy  on  deconffic- 
tion  of  aircraft  and  artillery  fire  is  a  prudent  course  of  action.  The  high  expected  cost 
of  flying  through  artillery  fires  must  be  matched  by  an  even  higher  expected  gain 
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in  order  for  the  mission  to  be  worth  the  possible  loss.  Furthermore,  given  current 
stand-off  bring  capabilities  and  high  altitude  ISR,  it  may  be  possible  to  adjust  flying 
missions  to  achieve  deconfliction  with  minimal  degradation  of  the  desired  mission 
effect.  Decision  analysis  can  aid  planners  in  trade-off  analysis  for  specific  scenarios 
given  the  probability  calculations  presented  in  this  paper,  but  in  general,  this  model 
indicates  deconfliction  is  the  recommended  course  of  action. 

3. 6  Conclusion 

Despite  the  many  collision  models  presented  in  the  literature,  the  gas  model  is 
the  best  probabilistic  model  for  addressing  the  problem  of  aircraft  fratricides  from 
artillery  fire.  It  provides  good  fidelity  with  low  data  requirements.  The  adjust¬ 
ments  that  are  necessary  to  apply  the  gas  model  to  the  fratricide  problem  require 
mathematical  contributions  from  diverse  fields.  Adjusting  the  “minefield”  model  to 
account  for  three  dimensions,  non-constant  shell  velocity  that  follows  ballistic  equa¬ 
tions,  and  integrating  for  uncertainty  in  the  input  parameters  all  combine  to  form 
a  unique  model  that  provides  an  original  contribution  to  the  literature  of  collision 
models. 

Numerical  results  from  the  model  indicate  unexpectedly  high  collision  prob¬ 
abilities  in  normal  scenarios.  Deconfliction  is  a  prudent  response  to  mitigate  the 
expected  costs  to  joint  campaigns,  but  decision  analysis  tools  may  be  necessary  in 
certain  situations  to  evaluate  the  best  course  of  action.  Sensitivity  analysis  on  the 
model  serves  the  dual  purpose  of  providing  validation  for  the  model  and  simultane¬ 
ously  suggesting  options  for  lowering  collision  risks. 

There  are  several  areas  where  this  collision  model  can  be  improved.  First, 
the  unexpectedly  high  collision  probabilities  suggest  that  a  simulation  model  could 
reasonably  provide  further  validation  of  the  model  without  running  the  risk  of  poor 
estimates  resulting  from  rare  occurances  due  to  tail  probabilities.  Better  ballis- 
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tic  models  that  incorporate  drag  coefficients  may  provide  better  estimates,  but  the 
marginal  improvements  in  accuracy  may  be  outweighed  by  the  increase  in  complex¬ 
ity.  Finally,  the  model  could  be  expanded  to  include  location  on  the  artillery  unit, 
target,  and/or  aircraft  position  using  inhomogeneous  intensity  parameter  values. 
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4.  Conclusion 


4-1  Research  Summary 

This  thesis  developed  a  useful  extension  to  the  gas  particle  collision  model.  It 
demonstrated  the  feasibility  to  account  for  spatial  dependencies  in  the  object  veloc¬ 
ities.  In  addition,  the  extended  Poisson  SPP  allows  the  analyst  to  compute  collision 
probabilities  when  the  density  of  objects  in  a  bounded  space  may  be  over-  or  under- 
dispersed  with  respect  to  the  Poisson  distribution.  This  even  allows  for  binomial 
or  negative  binomial  SPPs.  Finally,  multi-dimensional  birth-death  processes  allow 
the  modeler  to  represent  objects  entering  and  exiting  the  space.  Each  one  of  these 
enhancements  presents  an  original  contribution  to  the  literature  of  encounter  models. 

The  application  of  the  research  to  the  problem  of  aircraft  fratricides  from  ar¬ 
tillery  fire  provides  an  excellent  forum  to  test  the  capabilities  of  the  model.  Valida¬ 
tion  via  sensitivity  analysis  suggests  the  model  is  accurate.  The  application  provides 
useful  insights  to  the  problem  domain  that  can  affect  military  policy  and  doctrinal 
decisions. 

4-2  Future  Work 

Recent  theoretical  development  for  encounter  models  is  limited.  For  dynamic 
situations,  the  literature  is  limited  to  three-dimensional  cases.  This  is  probably  due 
in  no  small  part  to  the  difficulty  in  evaluating  the  probability  density  function  for  the 
collision  angles  in  other  dimensions.  The  three  dimensional  case  can  be  solved  with 
simple  u  substitution,  as  shown  in  Chapter  3.  However,  there  are  many  applications 
that  could  benefit  from  an  encounter  estimate  for  a  dynamic  two-dimensional  model. 
Higher  dimensional  models  may  also  be  useful,  though  their  applications  may  be 
more  limited.  Of  particular  interest  is  the  case  in  dimensions  higher  than  eight.  This 
is  because  it  has  been  shown  that  the  surface  area  of  a  unit  hypersphere  is  maximal 
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in  seven  dimensions  (among  integer  dimensions)  (Weisstein,  2009).  For  dimensions 
greater  than  seven,  the  surface  area  monotonically  decreases.  This  may  have  some 
interesting  effects  on  the  probability  of  the  angles  between  different  motion  vectors 
according  to  the  development  in  Chapter  3. 

Another  theoretical  development  that  is  still  lacking  in  the  literature  of  en¬ 
counter  models  is  that  of  inhomogeneous  spatial  distributions.  Mathematically  this 
would  probably  entail  some  sort  of  Gaussian  distribution  around  the  expected  loca¬ 
tion  of  the  object  at  any  given  time  (Fraser,  1951).  Naturally,  this  must  also  take 
into  account  the  probability  density  function  of  motions  of  the  objects.  If  the  objects 
move  in  uniformly  distributed  directions,  any  inhomogeneities  in  the  SPP  density 
function  will  dissipate  as  the  system  evolves.  One  means  to  maintain  the  spatial 
inhomogeneities  would  be  to  model  the  point  arrivals  using  the  spatial  birth-death 
process  and  distribute  the  arrivals  using  a  Wiener  process  for  the  time  since  their 
arrivals.  Similar  techniques  may  also  assist  in  estimating  uncertainty  in  the  defined 
motion  paths  of  the  encounter- er. 

The  application  presented  in  Chapter  3  may  be  improved  by  allowing  inho¬ 
mogeneous  spatial  density  functions.  This  would  reflect  more  realistic  conditions 
where  there  is  a  higher  probability  of  shells  or  artillery  to  exist  in  certain  regions 
(i.e.,  along  the  gun-to-target  line,  or  in  a  particular  orbit).  As  stated  in  the  pre¬ 
vious  paragraph,  this  would  necessitate  further  consideration  on  the  assumption  of 
uniformly  distributed  directions  of  motion. 

The  insight  that  collision  risks  are  non-negligible  highlights  the  need  for  opti¬ 
mization  techniques  to  optimize  the  management  of  a  congested  airspace  to  allow  for 
maximum  mission  effectiveness  with  minimum  risk.  This  may  involve  optimization 
methods  to  deconflict  the  airspace.  Alternatively,  it  may  be  desirable  to  assign  mo¬ 
tion  paths  through  a  conflicted  airspace  that  travel  to  designated  destination  points 
with  minimal  risk  of  fratricide. 
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Finally,  the  model  described  here  can  be  applied  to  a  number  of  other  problem 
scenarios  other  than  the  one  described  in  Chapter  3.  One  application  would  be 
dynamic  spatial  birth-death  systems  that  represent  mobile  ad-hoc  wireless  networks. 
Since  wireless  signals  require  spatial  locality  to  maintain  sufficient  strength,  it  would 
be  of  interest  to  know  the  probability  of  a  particular  node  going  out  of  signal  range, 
or  simply  dropping  off  the  network.  The  implications  of  such  estimates  may  assist 
in  developing  robust  policies  for  network  information  flow  that  can  maximize  the 
probability  that  certain  packets  can  reach  certain  nodes. 
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Appendix  A.  Blue  Dart 

Foregoing  Deconfliction:  Consider  the  Risk! 

Capt  Timothy  Holzmann  and  Dr.  Jeffery  Cochran 
Dept,  of  Operational  Sciences 
Air  Force  Institute  of  Technology 

Integration  of  air -to- surface  fires  and  surface-to-surface  indirect  fires  requires  appli¬ 
cation  of  appropriate  restrictions:  altitude,  time,  or  lateral  separation  -  (JP  3-09.34) 

“From  Deconfliction  to  Integration  to  Interdependence!”  is  the  current  mantra. 
Current  operations  tempos  are  driving  commanders  to  make  ”  Go”  or  ”  No  Go”  deci¬ 
sions  in  very  short  amounts  of  time,  while  simultaneously  increasing  joint  operational 
effectiveness  by  having  all  service  branches  work  together  to  achieve  a  synergistic  ef¬ 
fect.  For  example,  a  commander  may  want  Army  artillery  and  Air  Force  close  air 
support  to  work  together  against  some  target.  This  naturally  may  present  a  risk  to 
the  aircraft  if  it  were  to  fly  in  the  path  of  an  artillery  shell.  In  the  military  world, 
taking  unnecessary  risks  can  lead  to  deadly  consequences. 

Joint  doctrine  requires  commanders  to  eliminate  the  risk  that  the  artillery 
will  accidentally  hit  the  aircraft  by  instructing  military  planners  use  deconflict  the 
airspace:  designate  separate  times/places  for  the  artillery  units  to  shoot  and  the 
aircraft  to  fly  (Joint  Publication  3-09.3).  However,  recent  military  operations  in 
Iraq,  Afghanistan  and  exercises  on  the  Korean  peninsula  have  highlighted  several 
drawbacks  of  requiring  deconfliction.  It  can  (1)  limit  the  commander’s  ability  to 
respond  in  a  timely  manner,  (2)  prevent  him  from  massing  his  forces  effectively,  (3) 
limit  the  pilot’s  tactical  options,  (4)  provide  sanctuary  to  the  enemy,  and  (5)  hinder 
the  prosecution  of  time  sensitive  targets.  The  commander  is  left  in  a  bind:  accept 
no  risk  by  deconflicting  and  accept  these  limitations,  or  avoid  the  limitations  but 
accept  an  unknown  risk  by  flying  our  aircraft  through  our  own  artillery  fire.  Should 
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a  commander  take  the  risk?  Can  he  count  on  the  “Big  Sky  -  Little  Missile”  theory  to 
come  through,  or  should  he  bet  on  “Murphy’s  Law”  that  something  will  go  wrong? 

In  2001,  the  National  Defense  Industrial  Association  (NDIA)  performed  a  study 
to  test  the  “Big  Sky  -  Little  Missile”  theory  using  a  series  of  simulations  (Integra¬ 
tion,  Control  and  Deconfliction  of  Joint  Fires  Study,  Volume  1).  Their  simulations 
of  4  aircraft  flying  through  a  barrage  of  artillery  fire  generated  only  1%  collision 
probabilities.  They  conclude:  “Chance  for  conflict  still  a  risk  though  probability 
low.” 

More  recently  we  teamed  up  with  the  Warfighter  Applications  team  from  the 
Air  Force  LeMay  Center  to  build  mathematical  formula  for  the  question  of  ground- 
to-air  fratricide  risk.  Simulations  can  vary  results  from  one  run  to  the  next,  whereas 
a  formula  never  changes  its  answer.  However,  mathematical  models  are  frequently 
much  harder  to  construct.  In  this  case,  the  calculations  involved  ranged  from  high 
school  physics  formulas  to  graduate  school  mathematics.  The  Air  Force  Studies  and 
Analyses  agency  provided  preliminary  review  and  validation  of  the  model.  After 
further  consultations  with  experienced  military  operators,  our  AFIT  team  generated 
several  realistic  scenarios  involving  close  air  support  operations  and  unmanned  aerial 
vehicle  surveillance  operations.  Their  results  indicated  that  collision  probabilities  for 
very  normal  situations  can  range  up  to  1%  to  5%.  While  risk  values  depend  greatly 
on  the  particular  scenario,  both  the  NDIA  and  AFIT  studies  show  that  they  are 
frequently  higher  than  what  might  be  assumed  under  the  “Big  Sky  -  Little  Missile” 
theory. 

Probabilities  must  be  multiplied  by  the  high  price  tags  of  modern  aircraft  to 
generate  the  expected  cost.  Flying  a  $12  million  aircraft,  such  as  the  A-10  Warthog, 
through  artillery  fire  at  1%  chance  of  it  getting  hit  is  like  throwing  away  $120,000 
-  and  that’s  not  counting  loss  of  public  confidence  or  a  pilot’s  life.  Furthermore, 
what  would  the  loss  of  that  A-10  mean  to  the  unit’s  effectiveness  until  it  could  be 
replaced?  As  one  former  A-10  pilot  stated,  “It  is  one  thing  to  be  crawling  around  in 
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the  mud  with  the  infantry  when  some  other  pilot  accidently  gets  shot  down  by  US 
artillery;  it  is  quite  another  to  be  the  pilot  who  has  to  walk  home”  (Capt.  Hensley, 
“A  Fly  Paper”). 

We  cannot  turn  these  decision  processes  into  just  a  game  of  numbers.  Finding 
a  numerical  risk  value  is  easy  compared  to  shouldering  the  responsibility  of  the 
consequences.  However,  we  can  take  a  cue  from  the  quantitative  estimates:  the 
numbers  suggest  that  the  current  policy  of  deconfliction  is  a  prudent  course  of  action. 

Every  day  letters  are  going  home  to  mothers,  fathers,  wives,  husbands,  and 
children  expressing  our  regret  on  the  part  of  a  grateful  nation  for  the  sacrifices  of 
brave  members  of  our  armed  forces.  The  hardest  ones  to  write  are  those  when  we 
killed  one  of  our  own. 
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Probabilistic  Estimation  of 
Rare  Random  Collisions  in  3-Si 


Appendix  B.  Poster  Chart 
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Appendix  C.  Excel  Model  User  Guide 

The  software  implementation  is  in  the  form  of  a  Microsoft®  Office  Excel®  spread¬ 
sheet  Microsoft®  Corporation  (2006).  There  are  three  separate  worksheets  within 
the  file  label  “Start  Here,”  ’’Inputs  &  Outputs”,  and  ’’GaussQuad”  respectively. 
Each  of  these  three  worksheets  is  treated  below. 

C.l  Worksheet  1:  “Start  Here” 

The  intent  of  the  software  tool  is  to  present  the  minimum  level  of  complexity 
required  for  an  analyst  to  perform  the  necessary  calculations  and  find  the  desired 
results.  This  first  worksheet  provides  the  most  simple  presentation.  The  worksheet 
is  divided  into  several  regions:  the  model  calculator,  the  converter,  the  data  sugges¬ 
tions,  and  the  charts. 

C.t.l  The  Model  Calculator 

The  calculator  is  the  main  area  of  user  interface.  There  is  an  input  section,  an 
output  section,  and  a  notes  section.  Each  of  these  is  explained  below.  Figure  C.l 
displays  these  three  sections  as  they  are  laid  out  in  the  software  tool.  The  units  for 
each  measure,  whether  input  or  output,  are  designated  beside  the  appropriate  cell. 
Cells  intended  for  user  input  are  colored  white,  cells  with  calculated  values  have  a 
light  shade,  and  the  other  cells  have  a  darker  shade. 

The  inputs  are  grouped  by  the  component  of  the  scenario  that  the  respective 
parameters  describe.  The  artillery  section  requires  the  muzzle  velocity,  angle  of  fire, 
and  the  number  of  rounds  per  minute  that  the  artillery  unit  will  be  firing.  As  a 
simplifying  feature,  the  model  eliminates  the  need  for  further  inputs  by  assuming 
that  the  standard  conditions  of  the  U.S.  Army  field  manual  6-40  applies  HQ  (1999). 
These  conditions  are  reproduced  in  Figure  C.2.  In  addition,  the  model  assumes  that 
there  is  no  drag  on  the  projectile,  so  the  shell  follows  an  ideal  ballistic  arc.  Aircraft 
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_ Inputs 

Artillery  Unit 


Muzzle  Velocity 


Maximum 

733 

ft/sec 

✓ 

Minimum 

732 

ft/sec 

</ 

Angle  of  Fire  (0  =  horz,  90  =  vert) 

Maximum 

41 

0 

✓ 

Minimum 

40 

0 

✓ 

Rounds  per  minute 

2 

rpm 

✓ 

Aircraft 

Speed  (ground  speed) 

Altitude  (Above  Ground  Level) 
Maximum 
Minimum 
Number  of  aircraft 


Scenario 


Surface  Area 

4 

nm! 

Min  desired  miss  distance 

100 

ft  «/ 

Vulnerable  time  window 

10 

min 

Outputs 

P(frat)  (per  aircraft) 

Collision  rate 

E(#  frats) 

Var(#  frats) 

Mean  ±  1  St.  Dev  (#  frats  during  time  window) 

N=  0  P(0  frats)= 

197/313 

0.4963045 

4.9630448 

4.9630448 
[0,  0.877] 
1/143 

(=  62.94%) 
frats/min 

frats 

frats 
(=  0.7%) 

Notes 

None 

40  ft/sec 

101  ft  «/ 

100  ft  ✓ 

5  «/ 


Figure  C.l  A  screen  shot  of  the  Excel  model 
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STANDARD  CONDITIONS 

1 

WEATHER 

2 

AIR  TEMPERATURE  100  PERCENT  (59°F) 

3 

AIR  DENSITY  100  PERCENT  (1,225  gm/rrO 

4 

NO  WIND 

POSITION 

1 

GUN,  TARGET.  AND  MDP  AT  SAME  ALTITUDE 

2 

ACCURATE  RANGF 

3 

NO  ROTATION  OF  THE  EARTH 

MATERIAL 

1 

STANDARD  WEAPON,  PROJECTILE,  AND  FUZE 

2 

PROPELLANT  TEMPERATURE  (70°F) 

3 

LEVEL  TRUNNIONS  AND  PRECISION  SETTINGS 

4 

FIRING  TABLE  MUZZLE  VELOCITY 

5 

NO  DRIFT 

LEGEND:  gm/m3  -  grams  per  cubic  meter 

Figure  7-1.  Standard  Conditions. 


Figure  C.2  Excerpt  from  FM  6-40  containing  standard  firing  conditions 

are  described  by  their  altitudes,  speed,  and  quantity.  The  scenario  specifications 
include  the  size  of  the  operational  area  (ie.  a  top-down  view  of  the  area  where  the 
aircraft  are  flying  and  artillery  are  bring),  duration  of  the  airspace  conhiction,  and 
the  minimal  desired  miss  distance.  All  inputs  have  a  “sanity  check”  beside  them 
which  is  green  if  the  parameter  is  within  allowable  constraints,  yellow  if  it  is  on  the 
border,  and  red  if  it  is  outside  allowable  constraints.  For  example  negative  values 
will  be  marked  red.  These  bags  are  also  highlighted  in  the  notes  section  below. 

The  inputs  for  angle  of  bre,  muzzle  velocity,  and  aircraft  altitude  are  required 
in  ranges.  The  actual  values  will  be  assumed  to  be  uniformly  distributed  over  the 
specibed  ranges.  Ranges  must  be  of  positive  length  (ie.  max  must  be  greater  than 
the  min)  or  else  the  calculator  will  provide  erroneous  results. 

The  outputs  of  the  software  model  (seen  in  the  middle  box  of  Figure  C.l) 
include  a  probabilistic  assessment  of  the  risk  of  one  or  more  collisions  during  the 
conhicted  time  period.  This  is  the  basic  computational  result  from  Chapter  3.  In 
addition,  the  model  provides  the  rate  at  which  collisions  are  expected  to  occur, 
the  expected  number  of  collisions  over  the  whole  scenario,  and  its  variance.  Since 
variance  is  difficult  to  interpret  to  a  operational  community,  a  range  estimate  is  also 
provided  that  gives  reasonable  upper  and  lower  bounds  on  the  expected  number  of 
collisions  during  the  given  time  period.  Finally,  a  binomial  distribution  is  used  to 
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test  the  probability  that  a  given  number  of  fratricides  (up  to  the  number  of  aircraft) 
will  occur. 

The  notes  section  flags  important  notes  to  the  analyst  on  the  whole  scenario. 
Specifically  the  calculations  are  made  from  considering  9  different  combinations  of 
possible  angles  of  fire  and  muzzle  velocities.  The  notes  section  flags  to  the  user  cases 
where  the  projectile  fails  to  reach  the  minimum  and/or  maximum  altitudes  at  which 
the  aircraft  are  flying.  If  it  fails  to  reach  the  minimum  altitude,  the  probability 
of  collision  is  zero.  If  it  fails  to  reach  the  maximum  altitude,  the  probability  may 
be  also  reduced  from  the  case  where  it  passes  through  the  full  range  going  up  and 
coming  back  down. 

C.  1.2  The  Unit  Converter 

Also  provided  for  the  user  are  a  unit  conversion  table  to  aid  the  analyst  in 
deciding  the  right  inputs.  This  converter,  shown  in  Figure  C.3,  is  displayed  just 
to  the  right  of  the  model  calculator.  The  first  two  sections  of  the  converter  will 
calculate  the  appropriate  muzzle  velocity  of  artillery  fire  for  a  given  angle  of  fire 
and  distance  to  target;  and  it  will  calculate  the  appropriate  angle  of  fire  for  a  given 
muzzle  velocity  and  distance  to  target.  Below  this  are  various  unit  conversions  with 
the  units  specified. 

C.1.3  The  Data  Suggestions 

The  default  data,  as  seen  in  Figure  C.4,  gives  the  user  some  reasonable  inputs 
for  aircraft  velocities  and  altitudes.  The  data  provided  is  unclassified,  gathered 
from  Jane’s  military  resources  Daly  (2008),  and  Wikipedia  Wikipedia  Contributors 
(2009a, f,e,b,c,d).  The  intention  is  to  make  the  model  simple  to  use  for  the  average 
analyst. 
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Useful  calculations 

Find  Artillery  Velocity 

Distance  to  target 

25km 

Angle  of  fire 

45degrees 

Muzzle  velocity 

1625.129ft/sec 

Find  Artillery  Angle 

Distance  to  target 

10km 

Muzzle  velocity 

1625.129ft/sec 

Angle  of  fire 

11.78909degrees 

Conversions 

m2—> nm2 

lmiles2  = 

1.324341nm2 

miles  —>ft 

0.060386miles  = 

318.8406ft 

ft  —>  miles 

4921ft 

0.932008miles 

knots  —>  ft/sec 

330knots  = 

556.9758ft/sec 

mph  ->  ft/sec 

75mph 

llOft/sec 

mils  —>  degrees 

1200mils  = 

67.5° 

meters  —>ft 

4000meters  = 

13123.2ft 

Figure  C.3  The  conversion  table 


Baseline  Aircraft  Data 

Speed 

mph) 

Type 

Aircraft 

Cruise  alt 

Cruising  Maximum 

A-10 

340 

518 

F-15 

900 

Tactical 

F-16 

1320 

F-22 

1043 

Cargos 

C-17 

403 

C130 

400 

28000 

RQ-7A  Shadow 

75 

141 

Silver  Fox 

81 

12000 

UAVs 

FQM-151A  Pointer 

22 

56 

12500 

AV  RQ-11  Raven 

40 

60 

100 

RQ-1  Predator 

84 

138 

25000 

Figure  C.4  Table  of  Baseline  Data 
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C.l-4  Graphs 

The  graphs  provide  the  analyst  with  an  understanding  of  how  the  risk  varies 
over  various  uncertainties  in  the  input  parameters.  All  three  graphs  plot  the  surface 
of  the  probability  of  no  collisions  against  angle  of  fire  and  muzzle  velocity.  They 
differ  in  the  altitudes  at  which  the  calculations  are  made.  Collectively  the  three 
graphs  show  the  surface  over  the  full  range  of  uncertainty  in  the  inputs. 

C.2  Inputs  &  Outputs 

The  second  worksheet  gathers  the  inputs  from  the  first  page  in  a  unified  form. 
It  performs  range  checks  on  the  inputs  to  ensure  that  they  meet  the  range  require¬ 
ments.  It  also  collects  the  results  data  from  the  third  worksheet  and  performs  basic 
calculations  to  prepare  them  for  output  to  the  first  worksheet. 

C.3  Gauss  Quad 

The  final  worksheet  performs  the  Gaussian  Quadrature  calculations  specified 
in  Chapter  3  of  this  thesis.  Table  C.l  outlines  the  role  of  each  column  in  computing 
the  numerical  integration.  Figure  C.5  gives  a  screenshot  of  the  worksheet. 


C-6 


Table  C.l  Clarification  on  column  meanings  for  “GaussQuad”  worksheet 


Columns  A-C: 
Columns  D-F: 

Column  G: 
Columns  H-J: 

Columns  K-N: 

Column  O: 

Columns  P-Q: 
Columns  R-S: 

Column  T: 

Column  U: 
Column  V: 

Column  W: 
Column  X: 

Column  Y: 
Columns  Z-AA: 


Tx  value  is  a  mapping  from  [x_min,  x_max]  to  [-1,1] 

Sx  value  is  a  rank  value  of  the  corresponding  tx  value. 

1  is  low,  2  is  mid,  3  is  high,  0  will  not  be  included 
in  integration 

Sample  #  value  (1-27).  0  values  will  not  be  included 
in  integration 

Actual  values  of  altitude  (y),  launch  angle  (theta), 

and  muzzle  velocity  (vm)  corresponding  to  the  tx  values  from 

Columns  A-C 

Weights  for  numerical  integration.  Total  weight  is  the  product 
of  the  other  three. 

Maximum  height  possible  for  the  given  muzzle  velocity 
and  launch  angle 

The  minimal  and  maximal  altitudes  where  collision  is  possible 
Calculations  to  determine  expected  relative  velocities  of  shell 
and  aircraft 

Duration  (in  seconds)  that  the  shell  spends  in  the  collision 
altitudes 

The  expected  density  of  artillery  shells  in  one  unit  of  airspace 
The  total  relative  volume  covered  by  aircraft  traveling  at  the 
relative  speed  for  the  duration  of  the  vulnerable  time. 

The  collision  intensity  parameter 

Probability  of  one  or  more  shells  in  a  single  aircraft’s  path  at  any 
point  in  time  during  the  vulnerable  time 
Weighted  probability  used  for  numerical  integration 
Flags  if  shell  does  not  attain  aircraft  min  (col  Z)  or  max 
(col  AA)  altitude. 
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Screenshot  of  “GaussQuad”  worksheet 


Figure  C.5 


& o dddddddddddddddddddddddddo 


c 


Bibliography 

Air,  Land,  Sea  Application  Center.  Kill  Box:  Multi-Service  Tactics,  Techniques,  and 
Procedures  for  Kill  Box  Employment.  Joint  Publication  3-09.24,  June  2005. 

Alexander,  B.  “Aircraft  Density  and  Midair  Collision,”  Proceedings  of  the  IEEE ,  58 
(3):377-381,  1970. 

Alfano,  S.  “Satellite  Collision  Probability  Enhancements,”  Journal  of  Guidance, 
Control,  and  Dynamics ,  0629(3)  :588  592,  2006. 

Alfriend,  K.  T.,  et  al.  “Probability  of  Collision  Error  Analysis,”  Space  Debris,  1: 
21-35,  1999. 

Baddeley,  A.,  et  al.,  editors.  Lecture  Notes  in  Statistics.  Springer  Science  +  Business 
Media  Inc.,  New  York,  2006. 

Ball,  F.  “A  note  on  variation  in  birth  processes,”  Mathematical  Scientist,  20(1): 
50-55,  1995. 

Ball,  F.  and  Lynn,  O.  “Optimal  vaccination  schemes  for  epidemics  among  a  popula¬ 
tion  of  households,  with  application  to  variola  rnior  in  Brazil,”  Statistical  Methods 
in  Medical  Research,  15:481-497,  2006. 

Ball,  F.  and  Donnelly,  P.  “Interparticle  Correlation  in  Death  Processes  with  Appli¬ 
cation  to  Variability  in  Compartmental  Models,”  Advances  in  Applied  Probability, 
19(4):755-766,  1987. 

Bethea,  D.  and  Herbranson,  T.  “Calculating  the  number  of  Aircraft  volume  boxes 
hit,  within  a  certain  airspace,  by  a  ballistic  following  a  simple  trajectory,” .  Personal 
correspondence,  25  July,  2008. 

Bethea,  D.  and  Herbranson,  T.  “Bullet  Background  Paper  on  Operations  Research 
Fratricide  Study,” .  Online  Posting  to  Air  Force  Knowledge  Now,  Warfighter  Ap¬ 
plications  COP,  11  2008. 


BIB-1 


Bird /Wildlife  Avian  Strike  Hazard  (BASH)  Team.  “United  States  Avian  Hazard 
Advisory  System,”,  2008.  Available  online  at:  usahas.com  (Accessed:  3/2/2009). 

Brix,  A.  and  Chadoeuf,  J.  “Spatio-temporal  modeling  of  weeds  and  shot-noise  G 
Cox  processes,”  Biometrical  Journal,  44:83-99,  2002. 

Brown,  T.  C.  and  Donnelly,  P.  “On  Conditional  Intensities  and  on  Interparticlc 
Correlation  in  Non-Linear  Death  Processes,”  Advances  in  Applied  Probability,  25: 
255-260,  1993. 

Brush,  S.  G.,  editor.  Kinetic  Theory,  Vol.  I:  the  Nature  of  Gases  and  Heat.  Pergamon 
Press,  Ltd,  Oxford,  1965. 

Burden,  R.  L.  and  Faires,  J.  D.  Numerical  Analysis,  6th  ed.  Thompson  Brooks/Cole, 
Belmont,  CA,  2005. 

Campbell,  M.  “Ball  hits  dove.  1-in- 13-million  odds.  Dove  not  amused.,”  The  Globe 
and  Mail,  F9,  31  March  2001. 

Carlton-Wippern,  K.  C.  “Analysis  of  Satellite  Collision  Probabilities  due  to  Tra¬ 
jectory  and  Uncertainties  in  the  Position/Momentum  Vectors,”  Journal  of  Space 
Power,  12(4):349-364,  1993. 

Carpenter,  J.  R.  “Non-Parametric  Collision  Probability  For  Low- Velocity  Encoun¬ 
ters,”  Advances  in  the  Astronautical  Sciences,  127(1) :227-242,  2007.  Argues  for 
Monte-Carlo  simulation  of  satellite  collisions. 

Chan,  K.  F.  “Collision  Probability  Analyses  for  Earth  Orbiting  Satellites,”  Advances 
in  the  Astronautical  Sciences,  96:1033-1048,  1997. 

Chan,  K.  “Short-term  vs.  Long-term  Spacecraft  Encounters,”  In  Proceedings  of 
the  AIAA/AAS  Astrodynamics  Specialist  Conference  and  Exhibit ,  Providence,  RI, 
August  2004. 

Chapman,  C.  R.  and  Morrison,  D.  “Impacts  on  the  Earth  by  Asteroids  and  Comets: 
Assessing  the  Hazard,”  Nature,  367:33-7,  1994. 


BIB-2 


Chiang,  Y.-J.,  et  al.  “Geometric  Algorithms  for  Conflict  Detection/Resolution  in 
Air  Traffic  Management,”  In  Proceedings  of  the  36th  Conference  on  Decision  & 
Control ,  San  Diego,  CA,  December  1997. 

Clausius,  R„  On  the  Mean  Lengths  of  Paths  Described  by  the  Separate  Molecules  of 
Gaseous  Bodies ,  135-147.  In  Brush  (1965),  1858. 

Comas,  C.  and  Mateu,  J.  “On  soft  and  hard  particle  motions  for  stochastic  marked 
point  processes,”  Journal  of  Statistical  Computation  and  Simulation ,  77:1091- 
1121,  2007.  doi:  10.1080/10629360600864217. 

Daly,  M.,  editor.  Jane’s  Unmanned  Aerial  Vehicles  and  Targets,  Issue  30.  Jane’s 
Information  Group,  Inc.,  Alexandria,  VA,  2008. 

DeAngelis,  A.  “Is  It  Really  Worth  Running  in  the  Rain?,”  European  Journal  of 
Physics,  8(3):201-2,  1987. 

Faddy,  M.  J.  “Non-linear  stochastic  compartmental  models,”  IMA  Journal  of  Math¬ 
ematics  Applied  in  Medicine  and  Biology,  2:287-297,  1985. 

Faddy,  M.  J.  “On  Variation  in  Birth  Processes,”  Advances  in  Applied  Probability, 
22(2):480-483,  1990. 

Faddy,  M.  J.  “On  Variation  in  Poisson  Processes,”  Mathematical  Scientist,  19(1): 
47-51,  1994. 

Faddy,  M.  J.  “Markov  process  modelling  and  analysis  of  discrete  data,”  Applied 
Stochastic  Models  and  Data  Analysis,  13:217-223,  1995. 

Faddy,  M.  J.  “Extended  Poisson  Process  Modelling  and  Analysis  of  Count  Data,” 
Biometrical  Journal,  39(4):431-440,  1997. 

Faddy,  M.  J.  “Markov  death  process  modelling  and  analysis  of  binary  data,”  Statis¬ 
tics  and  Probability  Letters,  40:9-13,  1998. 


BIB-3 


Faddy,  M.  J.  and  Smith,  D.  M.  “Extended  Poisson  process  modelling  of  dilution 
series  data,”  Journal  of  the  Royal  Statistical  Society,  Series  C:  Applied  Statistics , 
57:461-471,  2008. 

Faddy,  M.  J.  and  Smith,  D.  “Modeling  the  Dependence  between  the  Number  of 
Trials  and  the  Success  Probability  in  Binary  Trials,”  Biometrics,  61:1112-1114, 
2005. 

Faddy,  M.,  et  al.  “Likelihood-Based  Modeling  and  Analysis  of  Possum  Trapping 
Data,”  Journal  of  Agricultural,  Biological,  and  Environmental  Statistics ,  6:235- 
242,  2001. 

Fraser,  D.  A.  S.  “Generalized  Hit  Probabilities  with  a  Gaussian  Target,”  The  Annals 
of  Mathematical  Statistics,  22(2):248-255,  06  1951. 

Glen  Shilland,  M.  U.,  February  2009.  Personal  Interview. 

Gross,  D.  and  Harris,  C.  M.  Fundamentals  of  Queueing  Theory,  Third  ed.  John 
Wiley  &  Sons,  Inc.,  New  York,  1998. 

“Guidance,”  Scottish  Natural  Heritage,  n.  pag.,  2000.  Available  online  at:  http: 
/ /www. snh.org.uk/pdfs/strategy/reneable/COLLIS. pdf  (Accessed:  09/29/2008). 

Hagenaars,  T.  J.,  et  al.  “Spatial  heterogeneity  and  the  persistence  of  infectious 
diseases,”  Journal  of  Theoretical  Biology,  229:349-259,  2004. 

Hensley,  J.  L.  “A  Fly  Paper,”  Field  Artillery  Journal,  12-16,  May-June  1984. 
Available  online  at:  http://sill-www.army.mil/FAMAG/1984/MAY_JUN_1984/ 
MAY_JUN_1984_PAGES_12_16.pdf  (Accessed:  3/2/2009). 

HQ  Department  of  Army,  Linked  States  Marine  Corps.  Tactics,  Techniques,  and 
Procedures  for  Field  Artillery  Manual  Cannon  Gunnery  (with  Change  1).  Field 
Manual  6-40  Cl,  December  1999. 

Inselmann,  E.  H.  “The  Generalized  Zacks  Model,”.  Technical  Report  TP  12-77, 
1977. 


BIB-4 


Jenkin,  A.  B.  and  Gick,  R.  A.  “Collision  Risk  Posed  to  the  Global  Positioning 
System  by  Disposal  Orbit  Instability,”  Journal  of  Spacecraft  and  Rockets,  39:532 
538,  2002.  doi:  10.2514/1.15398. 

Johnson,  N.  L.,  et  al.  “NASA’s  New  Breakup  Model  of  EVOLVE  4.0,”  Advanced 
Space  Research,  28:1377-1384,  2001.  doi:  10.1016/S0273-1177(01)00423-9. 

Johnson,  R.  E.  Introduction  to  Atomic  and  Molecular  Collisions.  1939. 

Joint  Chiefs  of  Staff.  Joint  Tactics,  Techniques,  and  Procedurs  for  Close  Air  Support 
(CAS).  Joint  Publication  3-09.3,  September  2003. 

Kim,  C.  “The  Effect  of  Sensor  Performance  on  Save  Minefield  Tran¬ 
sit,”  Master’s  thesis,  Naval  Postgraduate  School,  2002.  Avail¬ 
able  online  at:  http://www.nps.edu/academics/gseas/usw/docs/ 

NPSResearchNewsFebruary2003UnderseaWarfareResearchCenter.pdf  (Accessed: 
3/2/2009). 

Loeb,  L.  B.  The  Kinetic  Theory  of  Gases.  McGraw-Hill  Book  Company,  Inc.,  New 
York,  1934. 

Maxwell,  J.  C.  Illustrations  of  the  Dynamical  Theory  of  Gases,  148-171.  In  Brush 
(1965),  1858. 

Microsoft®  Corporation.  Microsoft ®  Office  Excel ®  2007  (12.0.6331.5000)  SP1 
MSO  (12.0.6320.5000)  [computer  software].  Microsoft®  Corporation,  Redmond, 
WA,  2006.  (Microsoft  product  screenshot(s)  reprinted  with  permission  from  Mi¬ 
crosoft  Corporation.). 

Mpllcr,  J.  and  Waagepeterson,  R.  P.  Statistical  Inference  and  Simulation  for  Spatial 
Point  Processes.  CRC  Press  LLC,  Boca  Raton,  FL,  2004. 

Mpller,  J.  and  Waagepeterson,  R.  P.  “Modern  Statistics  for  Spatial  Point  Processes,” 
Scandanavian  Journal  of  Statistics,  34(4): 643-684,  2006.  Available  online  at:  http: 
/ /www. math. aau.dk/~jm/NordstatPaper. pdf  (Accessed:  3/2/2009). 


BIB-5 


Musselman,  D.  A.  “Joint  Fires  Battlespace  Deconfliction:  Doctrinal  Emphasis 
to  Eliminate  Airborne  Fratricide,”  Master’s  thesis,  Air  Command  and  Staff 
College,  Air  University,  2008.  Available  online  at:  https://www.afresearch.org/ 
skins/rims/q_mod_be0e99f3-fc56-4ccb-8dfe-670c0822al53/q_act_downloadpaper/ 
q_obj_l  3439698-  f81e-473d-a350-679e3eb59a8a/display.aspx?rs=enginespage 
(Accessed:  3/2/2009). 

Neuenswander,  D.  M.  “Position  Paper:  Risk  of  fratricide  with  fixed  wing  aircraft 
and  fires,”,  2008. 

Podlich,  H.  M.,  et  al.  “Likelihood  Computations  for  Extended  Poisson  Process 
Models,”  InterStat^  1999.  Available  online  at:  http://interstat.statjournals.net/ 
YEAR/1999/articles/9909001.pdf  (Accessed:  2/12/2009).  September  #1. 

Podlich,  H.  M.,  et  al.  “Semi-parametric  extended  Poisson  process  models  for  count 
data,”  Statistics  and  Computing ,  14:311-321,  2004. 

Pollack,  A.  F.,  et  al.  “Tomahawk  Deconfliction:  An  Exercise  in  System  Engineering,” 
John  Hopkins  Applied  Physics  Laboratory  Technical  Digest ,  18(l):79-89,  1997. 
Available  online  at:  http://www.jhuapl.edu/techdigest/tdl801/pollack.pdf  (Ac¬ 
cessed:  3/2/2009). 

Quigley,  J.  and  Walls,  L.  “Nonparametric  Bootstrapping  of  the  Reliability  Function 
for  Multiple  Copies  of  a  Repairable  Item  Modeled  by  a  Birth  Process,”  IEEE 
Transactions  on  Reliability ,  54(4),  2005. 

Reese,  S.  and  Sondow,  J.  “Parabola,”  MathWorld  -  A  Wolfram  Web  Re¬ 
source,  n.  pag.,  2009.  Available  online  at:  http://mathworld.wolfram.com/ 
UniversalParabolicConstant.html  (Accessed:  2/25/2009). 

Smith,  K.  A.  “Joint  Transformation  of  Aerial  Interdiction  by  Enhancing  Kill  Box 
Operations,”  Master’s  thesis,  Air  Command  and  Staff  College,  Air  LIniversity, 
2006.  Available  online  at:  https://research.maxwell.af.mil/papers/ay2006/acsc/ 
Smith.pdf  (Accessed:  3/2/2009). 


BIB-6 


Stoyan,  D.  “Hard  Problems  with  Random  Systems  of  Hard  Particles,”  2003.  Avail¬ 
able  online  at:  http://isi.cbs.nl/iamamember/CDl/abstracts/papers/3108.pdf 
(Accessed:  3/02/2009). 

Strike,  Land  Attack  and  Air  Defense  (SLAAD)  Committee  of  the  National  De¬ 
fense  Industrial  Association  (NDIA).  “Integration,  Control  and  Deconfliction 
of  Joint  Fires  Study,  Vol  1,”.  2001.  Available  online  at:  http://www.ndia. 

org/Divisions/Divisions/StrikeLandAttackAndAirDefense/Documents/Content/ 
Content  Groups /Divisionsl  /  Strike,  _Land_Attack_and_Air  .Defense /PDFsl5/ 
concept_of_ops.pdf  (Accessed:  2/25/2009). 

Toscas,  P.  J.  and  Faddy,  M.  J.  “Fitting  the  Extended  Poisson  Process  Model  to 
Grouped  Binary  Data,”  Computational  Statistics ,  20:595-609,  2005. 

Weisstein,  E.  W.  “Hypersphere,”  MathWorld  -  A  Wolfram  Web  Resource ,  2009.  n. 
pag.  (2/10/2009)  http://mathworld.wolfram.com/Hypersphere. 

Wikipedia  Contributors.  “A-10  Thunderbolt  II,”  Wikipedia,  The  Free  Encyclope¬ 
dia ,  2009a.  n.  pag.  (2/25/2009)  http://en.wikipedia.org/w/index.php?titlc=A- 
10_Thnnderbolt_II&oldid=273332097. 

Wikipedia  Contributors.  “C-130  Hercules,”  Wikipedia,  The  Free  Encyclope¬ 
dia ,  2009b.  n.  pag.  (2/25/2009)  http://en.wikipeclia.org/w/index.php?titlc=C- 
130  _Hercules&oldid=270629986. 

Wikipedia  Contributors.  “C-17  Globemaster  III,”  Wikipedia,  The  Free  Encyclo¬ 
pedia,  2009c.  n.  pag.  (2/25/2009)  http://en.wikipeclia.org/w/index.php?title=C- 
17_Globemaster_III&oldid=273224299. 

Wikipedia  Contributors.  “F-15  Eagle,”  Wikipedia,  The  Free  Encyclopedia,  n. 
pag.,  2009d.  Available  online  at:  http://en.wikipedia.Org/w/index. php?title= 
F-15\_Eagle\&oldid=273027584  (Accessed:  2/25/2009). 


BIB-7 


Wikipedia  Contributors.  “F-16  Fighting  Falcon,”  Wikipedia,  The  Free  Encyclopedia , 
n.pag.,  2009e.  Available  online  at:  http: / / en. wikipedia. org/w / index. php? title= 
F-16\_Fighting\_Falcon\&oldid=273246493  (Accessed:  2/25/2009). 

Wikipedia  Contributors.  “F-16  Raptor,”  Wikipedia,  The  Free  Encyclopedia,  n. 
pag.,  2009f.  Available  online  at:  http: / / en. wikipedia. org/w / index. php? title= 
F-22\_Raptor\&oldid=272969998  (Accessed:  2/25/2009). 


BIB-8 


